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Calculations of ground-state and excited-state properties of materials have been 
one of the major goals of condensed matter physics. Ground-state properties of 
solids have been extensively investigated for several decades within the standard 
density functional theory. Excited state properties, on the other hand, were rela- 
tively unexplored in ab initio calculations until a decade ago. The most suitable 
approach up to now for studying excited-state properties of extended systems is 
the Green function method. To calculate the Green function one requires the self- 
energy operator which is non-local and energy dependent. In this article we de- 
scribe the GW approximation which has turned out to be a fruitful approximation 
to the self-energy. The Green function theory, numerical methods for carrying out 
the self-energy calculations, simplified schemes, and applications to various systems 
are described. Self-consistency issue and new developments beyond the GW ap- 
proximation are also discussed as well as the success and shortcomings of the GW 
approximation. 
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I. INTRODUCTION 



The Hamiltonian for a many-electron system is given by (in atomic units Ti — m — e ~ 1) 






(1) 



where y is a local external potential such as the field from the nuclei. Solving the above Hamil- 
tonian has been a major problem in molecular and condensed matter physics. The Coulomb 
interaction in the last term makes the Hamiltonian difficult to solve. For a small system, such as 
an atom or a small molecule, it is possible to obtain the many-particle ground-state wavefunction 
by means of the configuration interaction (CI) method (Boys 1950, Dykstra 1988). In this method 
the wavefunction is expanded as a sum of Slater determinants (Slater 1930) whose orbitals and 
coefficients are determined by minimizing the total energy. A very accurate ground state wavefunc- 
tion and energy can be obtained. The computational effort, however, scales exponentially with the 
system size so that applications to large molecules or solids are not feasible. For excited states, the 
computational effort becomes very large already for a small system. Fortunately, in practice we are 
interested in quantities which do not require the full knowledge of the wavefunctions. For example. 



we are interested in the total energy, excitation spectra and expectation values of single-particle 
operators which can be obtained from the Green function described in a later section. 

Approximate theories are usually concerned with finding a good single-particle approximation 
for the Coulomb term. The earliest of these theories is the Hartree approximation (1928) where 
the non-local Coulomb term is replaced by an average local Coulomb potential (Hartree potential) 
from all the electrons. Although it gives reasonable results, due to a cancellation between exchange 
and correlation, the Hartree theory is not accurate enough in many cases. An extension of the 
Hartree theory which takes into account the fermionic nature of the electrons leads to the Hartree- 
Fock approximation (HFA) (Fock 1930) where in addition to the average local Coulomb potential 
there is a non-local exchange potential which reflects the Pauli exclusion principle. For systems 
with an energy gap in their excitation spectra, the HFA gives a qualitatively reasonable result. In 
fact the HFA works quite well for atoms but for insulating solids, the energy gap is in most cases 
too large. The reason for this can be traced back to the neglect of correlations or screening which 
is not too important in atoms but crucial in solids. The HFA already takes into account to a large 
extent correlation between electrons of the same spin since the Pauli exclusion principle (exchange) 
prevents them to get close together. Two electrons of opposite spin on the other hand are allowed 
to occupy the same single-particle state at the cost of a large Coulomb energy. Correlations keep 
electrons away from each other, creating a screening hole around each electron which reduces 
the interaction with the other electrons and thereby the Coulomb energy. The energy cost of 
transferring an electron from one site to a neighbouring site is substantially reduced by screening. 
In the tight-binding limit, i.e. for localized states, the energy gap is approximately given by U and 
this is essentially the effective Coulomb energy of the states that form the gap. Thus, correlation 
or screening reduces the gap from its Hartree-Fock value. In metals, the absence of correlations in 
the HFA leads to qualitatively wrong results such as zero density of states at the Fermi level due 
to a logarithmic singularity in the derivative of the single-particle spectra with respect to k at kp 
(see, e.g., Ashcroft and Mermin 1976). 

To simulate the effect of correlations. Slater introduced the Xa method where the exchange 
potential is modelled by a local potential of the form V^ — an}/^ ^ derived from the electron gas 
and scaled by a constant a to simulate correlations, n is the local electron density (Slater 1951a, 
1951b, 1953, 1974). This method has been quite successful in calculating ground-state properties 
and excitation spectra but it is semicmpirical. The Xa theory may be regarded as a precursor of 
the modern density functional theory (DFT) (Hohenbcrg and Kohn 1964, Kohn and Sham 1965) 
which has become a standard method for calculating ground-state properties of molecules and 
solids. Recent reviews on DFT may be found in Jones and Gunnarsson (1989) and Dreizler and 
Gross (1990). In DFT, the ground-state energy can be shown to be a functional of the ground-state 
density and satisfies the variational principle with respect to the density. The explicit form of the 
functional in terms of the density is not known and such an explicit functional may not exist. Using 
the variational property of the energy functional, one arrives at a set of single-particle equations, 
the Kohn-Sham (KS) equations (Kohn and Sham 1965), to be solved self-consistently: 

[_iv2 + T/^ + F-]0, =e,^, (2) 



= Ei^^ 



(3) 



where V^ and V^'^ are the Hartree and exchange-correlation potential respectively. In practical 
applications, the functional containing the effects of exchange and correlations is approximated by 
the local density approximation (LDA) where the density in the exchange-correlation potential of 
the electron gas is replaced by the local density of the real system (Kohn and Sham 1965). The KS 
eigenvalues £i have no clear physical meaning except for the highest occupied which corresponds 
to the ionization energy (Almbladh and von Barth 1985a). Although there is no theoretical justifi- 
cation, they are often interpreted as single-particle excitation energies corresponding to excitation 
spectra of the system upon a removal or addition of an electron. In many cases, in particular in sp 



systems, the agreement with photoemission spectra is quite good but as will be described in the 
next section, there are also serious discrepancies. 

A proper way of calculating single-particle excitation energies or quasiparticle energies (Landau 
1956, 1957) is provided by the Green function theory (Galitskii and Migdal 1958, Galitskii 1958). 
It can be shown that the quasiparticle energies Ei can be obtained from the quasiparticle equation: 



h^V2(r) + y«(r)]*,(r) + /dV S(r,r'; S,)*.(r') = ^.*.(r) 



(4) 



The non-local and energy dependent potential S, or the self-energy, contains the effects of exchange 
and correlations. It is in general complex with the imaginary part describing the damping of the 
quasiparticle. It can be seen that the different single-particle theories amount to approximating 
the self-energy operator E. Calculations of E are unfortunately very difhcult even for the electron 
gas. We must resort to approximations and this review describes the GW approximation (GWA) 
(Hedin 1965a) which is the simplest working approximation beyond the HFA that takes screening 
into account. 



A. Problems with the LDA 

The LDA has been very successful for describing ground-state properties such as total energies 
and structural properties. There is no clear theoretical justification why the LDA KS eigenvalues 
should give excitation energies, and even the exact V^"^ is not supposed to give the exact quasipar- 
ticle energies. Nevertheless, the LDA KS eigenvalues are often found to be in good agreement with 
the quasiparticle energies measured in photoemission experiments. Despite of its success, there are 
serious discrepancies already in sp-systems and they become worse in d- and f-systems: 

• The band width in Na is 10-20 % too large, 3.2 cV in LDA vs 2.65 eV experimentally (see, 
however, section V-A for the experimental band width). 

• The band gaps in semiconductors Si, GaAs, Ge, etc. are systematically underestimated, as 
much as 100 % in Ge. 

• The band width in Ni is ^ 30 % too large, 4.5 eV in LDA vs 3.3 eV experimentally (Hiifner 
et al 1972, Himpsel, Knapp, and Eastman 1979) 

• In f-systems, the LDA density of states is in strong disagreement with experiment. 

• In the Mott-Hubbard insulators of transition metal oxides the LDA band gap is much too 
small (Powell and Spicer 1970, Hiifner et al 1984, Sawatzky and Allen 1984) and in some 
cases the LDA gives qualitatively wrong results. For example, the Mott-Hubbard insulators 
CoO and the undoped parent compound of the high Tq material La2Cu04 are predicted to 
be metals (see e.g. Pickett 1989) 

• The magnetic moments in the transition metal oxides are systematically underestimated 
(Alperin 1962, Fender et al 1968, Cheetham and Hope 1983) 

• In alkali-metal clusters, the ionization energies calculated within the LDA are too low com- 
pared to experiment (Ishii, Ohnishi, and Sugano 1986, Saito and Cohen 1988). 

• In the LDA, the image potential seen by an electron in the vacuum far from a surface decays 
exponentially instead of the expected — l/4(z — zo) decay where z is the coordinate normal 
to the surface and zq is the position of the image plane (Lang and Kohn 1973). 

Strictly speaking, one should not blame the LDA for all of these discrepancies since many of 
them are related to excited state properties which are outside the domain of DFT. However, excited 
state properties of finite systems can be calculated within the time-dependent extension of DFT 
(Runge and Gross 1984, Gross, Dobson, and Petersilka 1996). 



B. Theories beyond the LDA 

When discussing theories beyond the LDA, a distinction should be made between theories which 
attempt to find better energy functionals but which he within DFT and those theories which 
attempt to mimic the self-energy in order to obtain better quasiparticle energies but which are 
then usually outside DFT. The GWA belongs to the latter. 



1. Gradient corrections 

Since the LDA is based on the homogeneous electron gas, it is natural to take into account 
the inhomogeneity in the charge density of real systems by including gradient corrections in the 
energy functional. There are two main approaches. The first is a semiempirical approach where 
the exchange-correlation functional is modelled by a functional containing parameters which are 
adjusted to give the best fit to the cohesive energies of a number of " standard molecules" . The most 
successful of these models are due to Becke (1988, 1992, 1996). The other approach attempts to 
calculate the coefficients in the gradient expansion from first principles (Langreth and Mehl 1983, 
Svendsen and von Barth 1996, Springer, Svendsen, and von Barth 1996, for a recent development 
see Perdew, Burke, and Ernzerhof 1997 and references therein). While in general gradient correc- 
tions give a significant improvement in the total energy (Causa and Zupan 1994, Philipsen and 
Baerends 1996, Dal Corso 1996) there is almost no major improvement for quasiparticle energies 
(Dufek et al 1994). 



2. LDA + U 

One of the problems with the LDA is the absence of orbital dependence in the exchange- 
correlation potential. Since the potential does not distinguish between orbitals with different 
m-quantum numbers, for systems containing a partially filled d- or f-shell one obtains a corre- 
sponding partially filled band with a metallic type electronic structure unless the exchange and 
crystal field splitting create a gap between the up and down channel. Thus, the late transition 
metal oxides, which arc insulators, are predicted to be metals by the (nonspin-polarizcd) LDA. In 
the LDA+U method (Anisimov, Zaanen, and Andersen 1991, Anisimov et al 1993, Lichtenstein, 
Zaanen, and Anisimov 1995) and its generalization (Solovyev, Dederichs, and Anisimov 1994, 
Solovyev, Hamada, and Terakura 1996), an orbital dependent potential U acting only on localized 
d- or f-states is introduced on top of the LDA potential. For Mott-Hubbard insulators or rare earth 
metal compounds where the partially filled 3d or 4f bands are split by the Coulomb interaction, 
forming the upper and lower Hubbard band, the LDA-I-U works reasonably well (Anisimov et al 
1997). The bandstructure, however, is unsatisfactory. Another problem with the method arises for 
systems with partially filled 3d shells which are metallic, like the transition metals. In this case, 
the LDA+U would produce unphysical results since it would split the partially filled band. 



3. Self-interaction correction (SIC) 

Apart from the problem with orbital dependence in the LDA, there is another problem associated 
with an unphysical interaction of an electron with itself. In DFT, only the highest occupied 
state is free from self-interaction but in LDA, there is in general self-interaction for all states. 
This self-interaction is explicitly subtracted out in the SIC formalism resulting in an orbital- 
dependent potential (Cowan 1967, Lindgren 1971, Zunger et al 1980, Perdew and Zunger 1981). 
Self-interaction is significant for localized states and it tends to zero for extended states since in 
the latter case the charge is spread over the crystal and therefore the Coulomb interaction of an 
electron with itself is of order 1/A^. Since self-interaction is usually positive, one would expect the 



LDA eigenvalues for localized states to be too high, as is indeed the case. For atoms, SIC therefore 
lowers the LDA eigenvalues giving better agreement with experiment. The orbital-dependent 
potential in SIC can describe Mott-Hubbard insulators although the bandstructure is not likely to 
be satisfactory. SIC predicts all the 3d monoxides to be insulators except VO which is correctly 
predicted to be a metal (Svane and Gunnarsson 1990, Szotek, Temmerman, and Winter 1993, Aral 
and Fujiwara 1995). SIC however fails to give a delocalized paramagnetic solution for the doped 
high Tc compounds and a similar problem is anticipated in LDA+U. For more itinerant systems, 
SIC does not give localized solutions, and it then reduces to the LDA. Accordingly, application to 
semiconductors, e.g. Si, does not increase the LDA gap. 



4- A generalized KS scheme 

A recent attempt to improve the LDA description of quasiparticle energies is to choose a non- 
interacting reference system having the same ground-state density as the real system, like in the 
conventional KS scheme, but with a non-local potential (Seidl et al 1996). Since the potential is 
non-local, the choice is not unique, different non-local potentials may generate the same ground- 
state density. A particular choice is a non-local screened exchange potential minus its local form. 
By definition, the functional from which this potential is derived, is zero at the correct density. The 
KS equations consist of the usual LDA exchange-correlation potential plus the chosen non-local 
potential. For semiconductors Si, GaAs, Ge, InP and InSb, this method improves the values of the 
band gap. Applications to other systems have not been performed so far. 

Another recent scheme proposed by Engel and Pickett (1996) incorporates part of the correlation 
energy into the kinetic energy functional which may be thought of as mass renormalization. This 
scheme is shown to improve the description of band gaps in silicon and germanium but it gives 
negligible correction for diamond and carbon. 



5. The optimized effective potential method 

A natural way of improving the LDA would be to calculate the exchange energy exactly and 
to generate a local exchange potential by taking a functional derivative of the exchange energy 
functional with respect to the density (Kotani 1995, Kotani and Akai 1996, Bylander and Kleinman 
1995a,b). The correlation energy can be approximated by the LDA value. The original idea of this 
method is due to Talman and Shadwick (1976) in their work on atoms where the Hartree-Fock 
total energy is minimized with orbitals restricted to be solutions to single-particle Hamiltonians 
with local potentials. The scheme is also known as the optimized effective potential method. 
Applications of this approach to several semiconductors and insulators (C, Si, Ge, MgO, CaO, 
and MnO) yields encouraging results regarding the band gaps which in most cases are improved 
from the corresponding LDA values. However, the scheme is still within the density functional 
formalism and it is intended to improve the energy functional rather than the self-energy. One 
advantage of this scheme is the possibility to systematically improve the energy functional. With 
regard to GW calculations, the scheme may provide better starting points than the LDA. The 
scheme may be extended to include correlations by using a screened interaction potential. 



C. Motivations for the GWA 

The theories described above have drawbacks when applied to calculating quasiparticle energies. 
Gradient corrections attempt to improve total energies but do not address the problem of improv- 
ing quasiparticle energies. Indeed, applications to transition metal oxides, where the gap is zero 
or grossly underestimated by the LDA, do not give any significant improvement. SIC theory only 
applies to localized occupied states and numerical calculations show that the resulting eigenvalues 



are too low. The LDA+U is designed for systems with localized states split by the Coulomb cor- 
relation, forming the upper and lower Hubbard band. Applications to transition metals, however, 
would lead to unreasonable results. The generalized KS scheme using screened exchange, like the 
other theories discussed above, has no energy dependence which can be important in some cases. 
Moreover, the choice of the non-local potential is rather arbitrary. Since the theory has not been 
applied extensively, it is difficult to judge its usefulness. The exact exchange approach should 
in principle improve total energies when correlations are also taken into account but it does not 
address the LDA problems with quasiparticles. 

The GWA is derived systematically from many-body perturbation theory. The form of the self- 
energy in the GWA is the same as in the HFA but the Coulomb interaction is dynamically screened, 
remedying the most serious deficiency of the HFA. The self-energy in the GWA is therefore non- 
local and energy dependent. The GWA is physically sound because it is qualitatively correct in 
some limiting cases (Hedin 1995) which allows applications to a large class of materials, metals or 
insulators. 



• 



• 



• 



• 



In atoms, screening is small and the GWA approaches the HFA which is known to work well 
for atoms. 

In the electron gas, screening is very important which is taken into account in the GWA and 
for semiconductors it can be shown that the GWA reduces the Hartree-Fock gaps. 

For a core electron excitation, the GWA corresponds to the classical Coulomb relaxation 
energy of the other electrons due to the adiabatic switching-on of the core hole potential, 
which is just what is to be expected physically. 

For a Rydberg electron in an atom, the GWA gives the classical Coulomb energy of the 
Rydberg electron due to the adiabatic switching-on of an induced dipole in the ion core. 

For the decay rate and the energy loss per unit time of a fast electron in an electron gas, the 
GWA gives the correct formula. 



The GWA has been applied with success to a wide class of systems ranging from simple metals to 
transition metals and their compounds. 

So far, the GWA has been applied mainly to calculate single-particle excitation spectra, but 
it is also possible to calculate the total energy (Lundqvist and Samathiyakanit 1969, von Barth 
and Holm 1997, Farid, Godby, and Needs 1990) and the expectation value of any single-particle 
operator in the ground state. 



D. A short historical survey 

The earliest attempt to include correlations beyond the HFA in the form of GW theory was 
probably the work of Quinn and Ferrell (1958) for the electron gas. Their calculations, however, 
were limited to states around the Fermi energy and several approximations were made. DuBois 
(1959a, 1959b) also calculated the self-energy of the electron gas within a GW type theory but his 
calculations were only for small values of the electron gas parameter r^ < 1 or for high densities 
since (47r/3)rf = p, where p is the electron density. His results have therefore received less attention 
because they are not directly relevant to real metals which have r^ ~ 2 — 5 (Al r^ ~ 2, Cs rg ~ 
5). The first full calculation of the self-energy within the GWA for the electron gas was performed 
by Hedin (1965a). He also showed in a systematic and rigorous way how the self-energy can be 
expanded in powers of the dynamically screened Coulomb interaction, with the GWA as the first 
term in this expansion. Later on, Lundqvist (1967a, 1967b, 1968) did extensive calculations of the 
self-energy of the electron gas for various densities and studied the spectral functions. Rice (1965) 
used a different version of what is conventionally known as the GWA, including vertex corrections 
(corrections beyond the GWA). His results are similar to those of Hedin. Later on, Mahan and 
his group (Mahan and Sernelius 1989, Frota and Mahan 1992) performed extensive self-energy 



calculations for the electron gas using various forms of the GW^-type approximations, studying the 
importance of vertex corrections. 

Due to computational difficulties, the GWA was not applied to real materials until the mid 
eighties starting with the work of Hybertsen and Louie (1985a, b, 1986, 1987b) on semiconductors 
with encouraging results. At about the same time, Godby, Sham and Schliiter (1986, 1987a, b, 1988) 
did the same calculations and their results are in good agreement with those of Hybertsen and 
Louie (1986). We should also mention an earlier calculation for diamond using the tight-binding 
approach by Strinati, Mattausch, and Hanke (1982) although it is superseded by later calculations. 
The good results for semiconductors encouraged further applications to more complicated systems, 
transition metals and their compounds, jellium surface, sodium clusters and to f-systems. By now 
the GWA has become a standard method for including correlations beyond the HFA. 



II. THEORY 

In this section we describe a brief summary of the Green function theory and derive the GWA. 
More details on the Green function theory may be found in standard text books on many-body 
theory (e.g. Nozieres 1964, Fetter and Walecka 1971, Inkson 1984, Mahan 1990) and in the review 
article by Hedin and Lundqvist (1969). 



A. The Green function and the self-energy 

To study the electronic excitation spectrum of a solid, one performs photoemission experiments 
where photons with a certain energy to are used as projectiles to knock out electrons. By mea- 
suring the kinetic energy of the photoemitted electrons along a certain direction k and using the 
conservation of energy and momentum, the excitation spectrum -E(k) of the solid can be obtained: 

Lo = K.E. + E{k) (5) 

A photoemission experiment then measures the excitation spectrum of a solid with the presence of 
a hole (occupied density of states). An inverse photoemission experiment or BIS uses electrons as 
probes to measure the unoccupied density of states or the excitation spectrum with an additional 
electron. 

In the limit of large kinetic energy ("sudden approximation") (Hedin and Lundqvist 1969) of 
the photoemitted or probing electron, the spectrum is directly related to the one-particle Green 
function which is defined as 

iG{x,x')^{N\T[i^{x)tl;\x')]\N) (6) 

{N\^P{x)^P^x')\N) for t > t' (electron) , 

'{N\4>''{x')4>{x)\N) for t<t' (hole) ^ ' 

\N) is the exact N-electron ground state, iIj{x) is a field operator in the Heisenberg representa- 
tion which annihilates an electron at a; = (r,i), and T is the time-ordering operator which arises 
naturally from the time development operator defined later in equation (O) . The physical interpre- 
tation of the Green function is that for i' > i it is is the probability amplitude that a hole created 
at X will propagate to x' and for t > t' , the probability amplitude that an electron added at x' will 
propagate to x. Thus, the Green function describes the photoemission and inverse photoemission 
processes. 

From the Green function we can obtain 

• The expectation value of any single-particle operator in the ground state. 

• The ground-state energy. 



• The one-electron excitation spectrum. 



In this review, we are mainly interested in the excitation spectra. The first and the second 
properties have not been explored for real systems. 

From the Heisenberg equation of motion for the field operator 



^'-^^[M:^)M 



where the Hamiltonian is given by 



(firip^x)ho{x)ip{x) + - £r(fr'ip\x)ilj\x')v{r - r')ij;{x')ii){x) 



H 



we obtain the equation of motion for the Green function: 



d 



^di ^ ^^^^^ 



G{x,x')— / dx" M{x,x")G{x" ,x') ^ 5(x — x') 
where the mass operator (Hartree potential + self-energy) M is defined to be such that 
dxiM{x, xi)G{xi,x') — —i d'^riv (r — ri) 

x(iV|T Ut(ri,i)^(ri,t)V-(r,t)V-t(r',t')l 1^) 



(8) 



(9) 



(10) 



(11) 



ho is the kinetic energy operator plus a local external potential. The quantity on the right hand 
side is a special case of a two-particle Green function: 



G2(l,2,3,4)-(z)2(7V|r ^(l)^(3)^t(4)^t(2) |7V) 



(12) 



where 1 = xi = (ri, ii) etc. 

The self-energy may be evaluated in at least two ways, either by using Wick's theorem (Wick 
1950, Fetter and Walecka 1971) or by Schwinger's functional derivative method (Schwinger 1951, 
Martin and Schwinger 1959). We follow the latter. This is done by introducing a time varying 
field (j>{r,t) which is used as a mathematical tool for evaluating the self-energy and it will be set to 
zero once the self-energy is obtained. Working in the interaction (Dirac) picture we have 



\ijD(r,t))^U{t,toMD{r,to)) 
The time development operator U is given by 



U{t,to) =rexp 



-I I dT(j}{T) 
'to 



(13) 



(14) 



4){t) = I d^r(t>{Y,T)i}'^^{Y,T)i}D{v,T) 

The relationship between operators in the Heisenberg and Dirac representations is 

V-(r,i) = f7t(i,0)^i3(r,t)f7(i,0) 
The field operator V'd satisfies 



(15) 



(16) 



(17) 



so it is the same as the unperturbed (0 — 0) Heisenberg operator. The Green function can now be 
written as 



zG(l,2) 



{N°\T C/((X3, -00)^^0 (l)V'L (2) |A^°) 



(7VO|C/(oo,-oo)|7VO) 
By taking functional derivative of G with respect to (f) we get 

SG{1,2) 



<50(3) 



G(l,2)G(3,3+)-G2(l,2,3,3+) 



(18) 



(19) 



Using the above result for G2 in the definition of M in equation (^_1|) the term GG gives the Hartree 
potential V^ and we define 

(20) 



(21) 







T. = M-V". 






The equation of motion for the Green function becomes 








'4-M^) 


G{x,x')- j dx"^{x,x")G{x' 


,x') 


= 5{x - x') 


where 






Ho^ho + V" + ^ 






Using the identity 






"^ (G-^G) - G-i ^'^ 1 '^'^"^ G - 

6(f> 54) 6(/) 


5G 
5(j) 




and evaluating SG^^/Scj), where from equation (^) 










G-^^^|-^-H,-J:, 







we get 



E(l, 2) = i I d3d4G(l, 3+)T4^(l, 4)A(3, 2, 4) 
W is the screened Coulomb potential 

1^(1,2)= /d3e-i(l,3)w(3-2), 



where V is the sum of the Hartree and the external potential: 

V ^V" + 4 
A is the vertex function 

(5G-i(l,2) 



A(l,2,3) 



5V{i) 
= (5(1 - 2)(5(2 - 3) + 



'?S(1,2) 
5V{'i) 



<5(l-2)^(2-3)+ / d(4567)^|MG(4,6)G(7,5)A(6,7,3) 



(22) 



(23) 



(24) 



(25) 



(26) 



(27) 



(28) 



(29) 
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The second line is obtained from equation (|2j) and the last line by using the chain rule 6Y,/5V = 
{5Y, / 5G){5G / 5V) and by using the identity in equation ( |23| ) and the definition of A. 
Fourier transformation of equation ( |2l| ) gives (with now set to zero) 

[lo - Ho (r)] G (r, r', u) - j d^r"^{r, r", Lo)G{r",r\ lo) = S{r ~~ r') (30) 

If Go is the Green function corresponding to S = 0, then we have the Dyson equation 

G = Go + GoSG (31) 

The first term Go(l, 2) is a direct propagation from 1 to 2 without exchange-correlation interaction 
and S contains all possible exchange-correlation interactions with the system that an electron can 
have in its propagation from 1 to 2. 

In practical application, Go corresponds to Hq = ijHartroo ^ yxc -^^j^gj-g ^xc jg ^^qj^q local and 
energy-independent exchange-correlation potential, e.g. V^^j^^. In this case, the Dyson equation 
becomes 

G = Go + GoAEG (32) 

where AE = Y. -¥""=. 



B. The polarization and response function 

The response function is an important quantity in the evaluation of the self-energy. It is related 
to the inverse dielectric function e~^ as follows: 



£-1 = 


U V 

54) 




= 


1 + 


5p 


= 


1 + ' 


SpSV 

""SV 6(j) 


Ril 


,2) = 


6pil) 



(33) 

OV 0(p 

The response function is defined as 

(34) 

which gives the change in the charge density upon a change in the external field. We note that the 
above response function is a time-ordered one which is related to the physical (causal) response 
function R^ by (Fetter and Walecka 1971) 

ReR{uj) = ReR"iuj), Imi?(w)sgnw = Imi?^(w) (35) 

The polarization function is defined as 

which gives the change in the charge density upon a change in the total (external + induced) field. 
Noting that 

pil) = -iG{l,l+) (37) 

we can write 
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P(l, 2) = -« / d3d4G(l, 3)A(3, 4, 2)G(4, 1+) (38) 

In summary we have 

e-i = l + z;i? (39) 

e^l-vP (40) 

R = P + PvR (41) 

= w + vRv (42) 

C. The Hedin equations 

Summarizing the results in the previous sections, we arrive at the well-known set of coupled 
integral equations (Hedin 1965a, Hedin and Lundqvist 1969). From equations (|2^), (pl|), (|2g|), and 
(|2|) we have 

2(1, 2) = i / d{34)G{l, 3+)W^(l, 4)A(3, 2, 4) (43) 



G(l, 2) = Go(l, 2) + y d(34)Go(l, 3)2(3, 4)G(4, 2) (44) 

A(l, 2, 3) = 6(1 - 2)5{2 - 3) + J '^(4567)g|ii||G(4, 6)G(7, 5)A(6, 7, 3) (45) 



Wil, 2) = -^(1, 2) + / d(34)f (1, 3)P(3, 4)14^(4, 2) (46) 

where P is given in equation (p8|). Like G, A and VF satisfy Dyson-like equations. Starting 
from a given approximation for E the above set of equations can be used to generate higher 
order approximations. Although the equations are exact, a straightforward expansion for the self- 
energy in powers of the screened interaction may yield unphysical results such as negative spectral 
functions (Minnhagen 1974, Schindlmayr and Godby 1997). In fact, the expansion itself is only 
conditionally convergent due to the long-range nature of the Coulomb potential. So far there is 
no systematic way of choosing which diagrams to sum. The choice is usually dictated by physical 
intuition. 



D. Quasiparticles 

From the classical theory of Green functions the solution to equation (pm can be written in a 
spectral representation: 

G(r,r',.) = !: *■'-■"';''->' (47) 

"^-^ uj - Ei{uj) 

where '^i are solutions to the quasiparticle equation: 

ifo (r) *. (r, w) + / Sr I](r, r', c^)*,(r', w) = E,{Lo)^,(y, ^) (48) 
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In a crystal, the index i may be associated with the Bloch wave vector and band index. The 
eigenvalues Ei are in general complex and the quasiparticle wavefunctions are not in general or- 
thogonal because S is not Hermitian but both the real and imaginary part of E are symmetric. 
Suppose at some lo — uJi we find that tOi = Re£'i(a;i). If ImEi^LOi) is small, then the imaginary 
part of G is expected to have a peak at this energy (quasiparticle peak) with a life-time given by 
l/lmEiluji). It may happen that u — Rei?i(w) is zero or close to zero at some other energies and 
if the corresponding IinEi{Lu) are small, then we get satellites. For a non-interacting system, E is 
Hermitian and therefore Ei is real so that the quasiparticle has an infinite life-time. 

The spectral representation can also be obtained directly from the definition of G by inserting 
a complete set of {N ± l)-electron states in between the field operators and performing Fourier 
transformation, keeping in mind that the field operators are in the Heisenberg representation, i.e. 
tpit) = exp(illt)v3(0)exp(-illt). 

G^r.r'..)^^ ^""^'-''rl + /" <iJ-'^'-''rl (49) 



LU — LU' — iS .1 ^ LU — UJ' + IC 



The spectral function or density of states A is given by 

A{y,y',uj) = Im G'(r,r',Li;)sgn(a; — /^) 

TT 



Y, h{v)K[v')5[u ~^l + e{N~l, i)] (50) 

i 

+ Y.p*(r)p,{v')5[Lo ^ii~e{N + l, i)] (51) 



where 



/i,(r) = (Ar-l,z|^(r,0)|7V) (52) 

p,{Y)^{N + l,i\'^\Y,0)\N) (53) 

\N ± 1, z) is the ith eigenstate of the A^ ± 1 electrons with an excitation energy 

e{N ±l,i) ^ E{N ±l,i) - E{N ±1) (54) 

which is positive and E{N ± 1) is the ground-state energy of the A^ ± 1 electrons. The quantity /i 
is the chemical potential 

^i = E{N + l)-E{N) 
= E{N)-E{N-1) + 0{1/N) (55) 

The physical meaning of the poles of G is therefore the exact excitation energies of the A^ ± 1 
electrons. Since the poles of G in equations ( B^) and ( 49) must be the same, it follows that the 



real part of Ei(uji) are also the excitation energies of the A^± 1 electrons. For a very large system, 
the poles are usually so close together that it is meaningless to talk about the individual poles, and 
in an infinitely large system, the poles form a branch cut. In this case, it is more meaningful to 
interpret the excitation spectrum in terms of quasiparticles with energies Re£'i(wi) and life-times 
l/ImEi{uj, 



From equation (31), the spectral function A is schematically given by 



AiLo)^-y^\lmG,{uj)\ 

i 

Iv- |ImE,(w)| 



ReAE,(a;)|2 + |ImS,(w)|2 
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(56) 



where Gi is the matrix element of G in an eigenstate i/^i of the non- interacting system H^. A is 
usually peaked at each energy Ei — Ei + HeAT,i(Ei) (quasiparticlc peak) with a life-time given by 
l/|Im 'E,i(Ei)\ and renormalization factor (weight of the Lorentzian) 



^ _ dRe Aj:,{E,y ' 



duj 



< 1 (57) 



At some other energies cUp, the denominator may be small and A{ujp) could also show peaks or 
satellite structure which can be due to plasmon excitations or other collective phenomena. 

If we start with a single-particle Hamiltonian which is in some sense close to the true interacting 
Hamiltonian, the quasiparticles of the former Hamiltonian are just a set of i5-functions centred 
at the single-particle eigenvalues. If the interaction is switched on, typically the delta functions 
become broadened since the single-particle states can now decay to other excitations and lose some 
weight which might appear as collective excitations or satellite structures. The term quasiparticlc 
is arbitrary but we usually refer to quasiparticlc as an excitation originating from a single-particle 
state and to satellite as an excitation not contained in the approximate non-interacting system. 

The quasiparticlc energy can also be calculated to first order in Ei — Si as follows: 

S, =e, +ReAI],(£;,) 

.. .^ / ^ .„ ,dReAT.Jei) 

OUJ 

= e, + Z,Re AJ:,{e,) (58) 



E. The GW approximation 

The GWA may be regarded as a generalization of the Hartree-Fock approximation (HFA) but 
with a dynamically screened Coulomb interaction. The non-local exchange potential in the HFA 
is given by 

occ 

E-(r, r')^-J2 ^k„(r)Vk„(r')«(r - r') (59) 

kn 

In the Green function theory, the exchange potential is written as 

I]^(r, r', t~t')^ iG{r, r', t - t')v{r - r')6{t - t') 

which when Fourier transformed yields equation ( p9| ) . The GWA corresponds to replacing the bare 
Coulomb interaction w by a screened interaction W : 

S(l,2) = zG(l,2)T4^(l,2) (60) 

This is physically well motivated especially in metals where the HFA leads to unphysical results 
such as zero density of states at the Fermi level, due to the lack of screening. Formally, the GWA 
is obtained by neglecting the second term in the vertex function in equation (p9|), i.e. setting 
A(l, 2, 3) = (5(1 — 2)(5(2 — 3). Fourier transforming equation ( |60| ) we get 

E(r, r', ^) = ^j duj'Gir, r', lj + ij')W{r, r', lj') (61) 

For a non-interacting Go the imaginary part of the correlation part of the self-energy can be 
evaluated explicitly: 

occ 

Im E^(r, r',uj<fi) = nJ2 '/'k„(r)VL(r') Im W'ir, r', ek„ - w)0(£kn - uj) (62) 



kn 
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IniE^(r,r',w>/^) = -^ J2 ^knir)iJ;l„ir')lmW%r,r',LO - ekn)e{u; - e^n) (63) 

kn 

where 

]¥" = ]¥ -V (64) 

is the frequency-dependent part of W. The above results are obtained by expressing G and W in 
their spectral representations. The spectral representation of G is given in equation (p9). For W 
it is given by 



J-oo UJ-Uj'-td Jo 



dJ ^ \ ' .; 65 



D is proportional to the imaginary part of W and defined to be anti-symmetric in w: 

D(r, r', iS) = -ilm VK'=(r, r', tj)sgn(w) (66) 

TT 

D(r,r',-c^) = -i?(r,r',c^) (67) 

The spectral representation of the correlation part of the self-energy is 

J-oo UJ - UJ' - 10 J^ LU-Uj'+ld 

where 

r(r,r',tj) = --ImE'=(r,r',tj)sgn(w-^) (69) 

TT 

The real part of S'^ can be obtained by performing the principal value integration (Kramers-Kronig 
relation or Hilbert transform). 

A physically appealing way of expressing the self-energy is by dividing it into a screened-exchange 
term Ssex and a Coulomb-hole term Scoh (COHSEX) (Hedin 1965a, Hedin and Lundqvist 1969). 
It is straightforward to verify that the real part of the self-energy can be written as 

occ 

RcEsEX (r, r', t^) = - ^ 4^^{vm {r')ReW (r, r\u ~ £,) (70) 



ReI]coH(r,r',c^)= Vv^(r)^*(r')P / 

"'0 



,^,Dir^rWl (,^) 



The physical interpretation of Scoh becomes clear in the static approximation due to Hedin 
(1965a). If we are interested in a state with energy E close to the Fermi level, the matrix element 
(V'IReEcoH(^)l^) picks up most of its weight from states with energies Si close to E in energy. 
We may then assume that E — Si is small compared to the main excitation energy of D, which is 
at the plasmon energy. If we set E — Si = 0, we get 



ReEcoH (r, r') = ^<5 (r - r') W^ (r, r'. 



0) (72) 



This is simply the interaction energy of the quasiparticle with the induced potential due to the 
screening of the electrons around the quasiparticle. The factor of 1/2 arises from the adiabatic 
growth of the interaction. In this static COHSEX approximation, Ecoh becomes local. 

The polarization function needed to evaluate W is calculated within the random phase ap- 
proximation (RPA) (Pines and Bohm 1952, Bohm and Pines 1953, Lindhard 1954, Pines 1961, 
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Gcll-Mann and Brueckner 1957) which corresponds to neglecting the second term in the vertex 
function and using a non-interacting Go : 



occ unocc 



spin kn k'n' 

'- -. '- ^1 (^3) 

We have used the fact that for every ip^n there is V'lkn with the same eigenvalue due to the time- 
reversal symmetry. The wavefunction ^kn is normahzed to unity in the entire space. The physical 
meaning of the RPA is that the electrons respond to the total field (external -I- induced field) as if 
they were non-interacting. 



III. NUMERICAL METHODS 

One of the main computational problems is to calculate the polarization in equation d73). Since 
P can be long-ranged for crystals, the conventional method of calculating it is to use Bloch basis 
functions. Noting that P(r -I- T, r'-f T) = P(r, r') , P can be expanded in general as follows 

P{ry,Lu) = ^i?q,(r)P,,(q,c^)i?;^.(r') (74) 

where {Bqi} is a set of Bloch basis functions large enough to describe P and the sum over q is 
restricted to the first Brillouin zone. The matrix elements of P in the basis are given by 

occ unocc 



Pij{q,Uj) = E E E (^q»V'kn|V'k+qn')(V^k+qn'|V'knB^ 
spin,k n n' 

1 1 



^ - £k+qn' + £kn + iS UJ + Sk+qn' - £kn - iS 

Similarly, the Coulomb potential v can be expanded as in equation (u4). The screened potential 
W can then be calculated from equation (|2|). 

The matrix element of the imaginary part of the self-energy in a state ipqn is given by 






ImE^„(w) = -^2^ 2^ 2^(^q«V'k-q,«'|-Ski) lmW,'j(k,uj- ek-q,„') 

k n'<fi ij 

X (Bkj|l/'k-q,n''0qn) 6'(^ -£k-q,n') ior U < fl (75) 

If we are only interested in quasiparticle energies, it is more favourable to perform the frequency 
integration in the expression for the self-energy in equation mU) along the imaginary axis (Godby, 
Schliiter, and Sham 1988). In this case, W must also be calculated along the imaginary axis which 
is advantageous since the pole structure along the real axis is avoided. For states away from the 
Fermi level, there is in addition a contribution to the self-energy from the poles of the Green 
function. 

The choice of basis functions depends on the type of materials we are interested in. For sp 
systems, plane waves are appropriate especially when used in conjunction with pseudopotentials. 
For systems containing rather localized states such as the 3d and 4f systems, a large number of 
plane waves would be needed and therefore localized basis functions are more suitable. 
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A. Plane-wave basis 

In this case i?kj -^ exp [i(k + G) • r] /v" where the index j is represented by the reciprocal 
lattice vector G and il is the unit cell volume. This is probably the simplest basis with the 
following advantages: 

• Programming ease: The matrix elements can be calculated easily particularly when the 
wave-functions are also expanded in plane waves. 

• The Coulomb potential is diagonal with matrix elements given by 47r/|k + G| . 

• Good control over convergence. 
The disadvantages are: 

• It is not feasible to do all-electron calculations. In many cases, it is essential to include core 
electrons. For example, the exchange of a 3d valence state with the 3s and 4p core states in 
the late 3d transition metals is overestimated by the LDA by as much as 1 eV which would 
lead to an error of the same order in the pseudopotential method. 

• The size of the response matrix becomes prohibitively large for narrow-band systems due to 
a large number of plane waves. 

• No direct physical interpretation. 



B. Localized basis 

For systems containing 3d or 4f electrons, plane- wave basis becomes very costly. Methods based 
on the linear-muffin-tin-orbital (LMTO) basis or Gaussian basis are more appropriate. In the 
LMTO method (Andersen 1975), the wave functions are expanded as follows 

V'k^W = ^XRL(r,k)6k„(RL) (76) 

RL 

where x is the LMTO basis which in the atomic sphere approximation (ASA) for r in the central 
cell is given by 

XRL(r,k) = 0RL(r) + Y, ^R'L' (i-)^R'L',RL(k) (77) 

R'L' 

</'RL(r) = 'i^'Riir)YL{Q) is a solution to the Schrodinger equation inside a sphere centred on an 
atom at site R for a certain energy e^/, normally chosen at the centre of the band, ^j^^, is the 
energy derivative of (Jjul at e^ . One advantage of the LMTO method is that 0rl is independent 
of k. When forming the polarization function in equation ( [73[ ) we have products of wave functions 
for r or r' of the form 

^^ :^ [(j)(j) + {(j} (I) + (I) (l))h+ (j)(l) h^]b^ (78) 

It is then clear that the sets of products (jjcj), (j) (/), (j)(j) form a complete basis for the polarization 
function and the response function (Aryasetiawan and Gunnarsson 1994a) since the latter can be 
written as 

R = P + PvP + PvPvP + ... (79) 

so that the basis for R is completely determined by that of P. Although this product basis is not 
complete for the Coulomb potential, it is of no consequence since v is always sandwiched between 
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two P's. It can be easily shown that the product basis is also complete for the self-energy. Thus 
schematically 



E(kn,co) = (V'knKGiy|^k„) 

= (V'knV'klV'V'kn) + i'iplcnM'^Rvl'Iplplcn) 



(80) 



which shows that it is sufficient to expand v in the product basis and the latter is therefore a 
complete basis for the self-energy. 




50,0 
u(eV) 



100.0 



FIG. 1. The energy loss spectra of Ni for q — (0.25, 0) 2n/a, a — 6.654ao. The large dots correspond 
to the experimental spectrum taken from Feldcamp, Stearns, and Shinozaki (1979). The solid line and small 
dots are respectively the loss spectra with and without local field corrections due to the inhomogeneity in 
the charge density. Both spectra are calculated with 4s, 4p, 3d, 4d, 4f, and 5g orbitals, including an empty 
sphere at (0.5, 0.5, 0.5)a and core excitations. After Aryasetiawan and Gunnarsson (1994). 

The number of product functions is still large, with nine spd orbitals we have 2[9(9 + l)/2] 
products of (j)(f> and (j)(j) and 9x9 products of (j)(f> giving in total 161 product functions. With spdf 
orbitals the number of product functions is 528. It can be reduced considerably without too much 

loss of accuracy by neglecting the (j) terms since they are small. Moreover, in the polarization 
function there are no products between conduction states. Therefore in sp and d systems products 
of 4'd4'd and (j)f(t>f respectively can be neglected. After these eliminations, the remaining product 
functions turn out to have a significant number of linear dependencies, typically 30-50 %, which 
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can be eliminated further giving ~ 100 optimized product functions per atom for spdf orbitals. A 
product function is given by 

Ba{r) = (/)RL(r)(/)RL/(r) 

= ^m (r)^Rr (r)Fi iQ)YL' (n) (81) 

where a = (R,LL'). This function is non-zero only inside an atomic sphere centred on atom R. 
There are no products between orbitals centred on different spheres. The optimization is performed 
by calculating the eigenvalues of the overlap matrix Oap — (BalBp) and subsequently neglecting 
eigenvectors with eigenvalues less than a certain tolerance. The optimized product functions are 
then linear combinations 

Bqi = y BjZ'fa (82) 

7 

where z are the eigenvectors of the overlap matrix: Oz — ez. Due to the localized property of 
the basis, the calculation of the dielectric matrix scales as N'^. This approach has been used with 
success to calculate loss spectra (Aryasetiawan and Gunnarsson 1994b) and self-energy for 3d 
systems (Aryasetiawan and Gunnarsson 1995). The loss spectra of Ni calculated using the product 
basis is shown in figure |l|. 

Another approach using a localized basis set is based on Gaussian functions (Rohlfing, Kriiger, 
and PoUmann 1993). The wave functions and the dielectric matrix are expanded in this basis. 
For Si, for example, the number of Gaussian orbitals needed to perform GFK calculation is 40-60 
whereas 350 plane waves are needed in the conventional approach. This approach has been applied 
with success to a number of semiconductors and insulators and to semicore states in Si, Ge, and 
CdS as well as to the Si surface. 



C. The plasmon-pole approximation 

One of the major computational efforts in self-energy calculations is the calculation of the 
screened interaction W. The physical features of W are well known; the imaginary part of W 
is characterized by a strong peak corresponding to a plasmon excitation at the plasmon frequency. 
This is particularly evident in the case of the electron gas or the alkalis such as Na and Al. The 
plasmon-pole approximation assumes that all the weight in Im W resides in the plasmon excita- 
tion (Lundqvist 1967, Overhauser 1971, Hedin and Lundqvist 1969, Hybertsen and Louie 1986). 
In the case of the electron gas, this is strictly true in the limit of long wave length q ^ 0. For 
finite q, the spectrum contains also particle-hole excitations at lower energies. The particle-hole 
spectrum eventually merges with the plasmon excitation as q gets larger. Thus, in the simplest 
form, the plasmon-pole approximation is given by (Lundqvist 1967, Hedin and Lundqvist 1969) 
Ime~^ (Q:^) — ^q^ {^ ^ ^q)- Thc two parameters Aq and ujq are determined from the static- and 
/-sum rules. In the generalized plasmon-pole approximation due to Hybertsen and Louie (1986), 
each matrix component of thc inverse dielectric function is written (for positive frequency) 

Im e^jQ, (q,w) == Agc (q) 6 {lj - ljgg' (q) ) (83) 

The corresponding real part is given by 



^Ic (q) 

^^ - ^GG' (q) 



Re e^\,, (q, Lo) = <5gg' + ^ _''''2 .^^ (84) 
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FIG. 2. Comparison between the numerically calculated inverse dielectric function (Walter and Cohen 
1972) and the corresponding plasmon-pole results for Si. After Hybertsen and Louie (1988a). 



The effective bare plasma frequency r^cG' is defined below. The unknown parameters Agg' (q) 
and cjgg' (q) are determined from the static limit of e^^: 



2 f°° 1 
Re e^^Q, (q, 0) = 6gg' + -P du;-lm e^\^, (q, 



u) 



and the /-sum rule: 



GG ^4' ^ 2 P |q+G|' p(0) 2 °° 



(85) 



(86) 



Thus, there are no adjustable parameters. The /-sum rule is true for the exact response function 
since it is obtained from the double commutator [[7J, pq+c] ,Pc^+G'] which yields 



\Y.^Es- Eo) { (0|/}q+Gk) {s\p[^^, |0) + c.c.} = (q + G) . (q + G') p (G - G') 



(87) 



The states \s) are the exact many-body excited states and p (G) is the Fourier component of the 
electron density. It can also be shown from the definition of the linear response dielectric function 
in terms of the ground-state matrix element of the commutator of the density operators that 
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f°° TT 

/ dtj oj Im Egg, (q,cj) = w (q + G) 

Jo 2 



:i ^ (ii;, - ii;o) {(0|/5q+Gk)(s|p]^+G,|0) + c.c.} 



The /-sum rule follows immediately from the last two equations (Hybertsen and Louie 1986) re- 
calling that ujp = 47rp(0). The quality of the plasmon-pole approximation is illustrated in figure 

i 

One drawback of the Hybertsen and Louie model is that the plasnion frequency may become 

complex, which is somewhat unphysical. A different model proposed by von der Linden and Horsch 
(1988) circumvents this problem by using the dielectric-bandstructure approach. Here one defines 
a Hermitian dielectric matrix 

47r 
|q + G| |q + G | 

«°GG^ (q, c.) ^ y y V T (M|e'('^^-)-1k + q,^)(k + q,^|e-('^^«')-1k,.) 
GG V4, y ^ /^ ^ ^ CJ - £k„ + ek+q,m ^ ' 

spin k n m ^' 

The inverse static dielectric matrix is expressed in its eigen representation 

oo 

eGG'(q.O) = -^GG' +Ec/q,:(G) [dr\cO - 1] C/;,(G') (91) 

The matrix U is formed by the eigenvectors of the inverse dielectric matrix. The plasmon-pole 
approximation is then obtained by introducing the frequency dependence in the eigenvalues: 

d-\ci,^)-l^ , "'i^l , (92) 

uj^ - [uji(q) - i6\ 

The eigenfunctions for lu ^ are approximated by the static eigenfunctions (von der Linden and 
Horsch 1988). As in the Hybertsen and Louie approach, the unknown quantities Zi and uji are 
determined from the static limit of e~^ and the /-sum rule. These give 

and 

-f (q) - 7^#rT (94) 

1-rf, (q) 

It can be shown that the pole strength Zi can be expressed in terms of the real-space eigenpotentials 

^.(q) = ^/rf'^p(r)|VvI/q,(r)p (95) 

where 

so that the pole strength is positive definite. The eigenvalues d~^ of the static dielectric matrix 
lies between (0,1) which implies that the plasmon frequencies cui are all real (von der Linden and 
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Horsch 1988). A generalization of the plasmon-pole model has been proposed by Engel et al (1991) 
and Engel and Farid (1992, 1993). 

One drawback of the plasmon-pole approximation is that the imaginary part of the self-energy 
is zero except at the plasmon-poles. As a consequence, the life-time of the quasiparticles cannot 
be calculated. Another drawback is its limited applicability to other than sp systems. For more 
complex systems where the plasmon excitations merge with the particle-hole excitations, it is not 
clear anymore if the plasmon-pole approximation is appropriate. 



D. The space-time method 

Conventional ways of performing self-energy calculations express all quantities in frequency and 
reciprocal space. This is natural because in solids the Bloch momentum is a good quantum number 
and by working in reciprocal space, one takes advantage of the lattice translational symmetry. In 
the extreme case of the electron gas, the self-energy becomes diagonal in k. In fact, it has been 
found numerically that in the Bloch-state representation, the self-energy is almost diagonal even 
in systems that bear little resemblance to the electron gas (Hybertsen and Louie 1986, Aryase- 
tiawan 1992a). The reason for this has been clarified recently by Hedin (1995). The frequency 
representation is also natural because most experiments generally focus on energy-dependent mea- 
surements. For example, photoemission spectra are measured as functions of momentum and 
energy. Theoretically, however, the self-energy becomes a multidimensional convolution when ex- 
pressed in momentum-energy space which is computationally very expensive as may be seen in 
equation (|75|). In space-time representation, on the other hand, the self-energy in the GWA takes 
a simple multiplicative form 

I](r,r',i)=zG(r,r',t)I^(r,r',t) (97) 

The response function also has a simple form 

P" (r, r', t) = -iG° (r, r', t) G° (r', r, -t) (98) 

and the non-interacting Green function G" is given by 

f *EknV'k„(r)^^„(r')e--*, t<0, 
G°(r,r',i) = <^ (99) 

l-*Ekr>k„(r)^L(r')e— *, t>0 

Evaluation of G° in real time is not advantageous due to the oscillatory exponential term. In 
imaginary time this exponential term decays rapidly and the summation can be performed easily. 
G" becomes (Rojas, Godby, and Needs 1995) 

f *Er„>k„(r)^L(r')e-^-^ r < 0, 
G" {r,r', ir) = { (100) 

l-*Ekr'^kn(r)^L(r')e-^-^ r>0 

This is obtained by analytically continuing the expression for G° in equation (^) from real to 
imaginary energy and taking Fourier transform from imaginary energy to imaginary time. From 
translational and lattice symmetry, r can be restricted to lie in the irreducible wedge but r' must 
run over all grid points which are located inside a large sphere of radius i?max centred on r. The 
method was tested for jelhum and Si and typical values for the grids are Ar = 0.5 a.u., At = 0.3 
a.u., i?max — 18 a.u., and T,„ax — 10 a.u., which give convergence in the quasiparticle energy 
differences to 0.05 eV and absolute values to 0.1 eV (Rojas, Godby, and Needs 1995). 

Although the polarization function P*^ takes a particularly simple form in space-time represen- 
tation, calculating e^^ in real space is still computationally prohibitive because of the large size of 
matrices in inverting the dielectric function. In practice, one performs six-dimensional fast Fourier 
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transform from (r,r') to (k, G, G') and one-dimensional Fourier transform from imaginary time 
to imaginary energy and solves for e~^ as a matrix equation in G, G' for each k and iu). After 
forming M^gg' (k, iuj) one Fourier transforms back to real space and to imaginary time representa- 
tion. The advantages of this method are that the computational effort is much reduced compared 
to the conventional techniques since the double summation over k points and bands are avoided 
and there is no need for a plasmon-pole approximation. 

For comparison with experiment, the self-energy has to be calculated at real frequencies. This 
is achieved by first calculating the self-energy and the matrix element of the self-energy correction 
directly in real space ("010115^ (*''") ~ ^™|V'kn) and then Fourier transforming the result from imag- 
inary time to imaginary energy. The imaginary energy representation can then be analytically 
continued to real energies by fitting the self-energy to the multipole form 



ao-f > 



(101) 



i=l 



where the parameters at and bi are in general complex. A good fit is obtained with only n = 2 
with an RMS error of 0.2 %. An example is shown in Fig. a. The advantage of calculating the 
self-energy along the imaginary axis is that one avoids the sharp pole structures in both G and W. 
For quantities which require frequency integration, it is then useful to perform the integration along 
the imaginary axis whenever possible. However, when the detailed structure of the self-energy is 
required along the real axis, knowledge of it at a few points along the imaginary axis is not likely 
to be sufficient. 
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FIG. 3. The real part of the matrix elements of the self-energy operator of Si continued onto the real 
ajds for the first 8 bands at k = 0. Inset: The matrix elements calculated along the imaginary axis for 
band 4 (the valence band maximum) together with the fitted form (with two poles). After Rojas, Godby, 
and Needs (1995). 

One important aspect of the space-time method is its scaling with respect to the system size (Ro- 
jas, Godby, and Needs 1995). It is found numerically that the range of W and E is rather material 
independent so that the parameters i?max and Ar do not change with system size. This means that 
the information stored as well as the computational time for the non-interacting response function 
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should scale linearly with the system size. However, the calculation of the dielectric function in- 
volves matrix inversions which scale as N^. It is still favourable compared with the conventional 
plane-wave basis approach which scales as N^ but comparable to the localized basis approach which 
also scales as N'^ (Aryasetiawan and Gunnarsson 1994b). 

An interesting mixed-space approach for calculating the polarization function was recently pro- 
posed by Blase et al (1995). The polarization function is written as 



P(r, y', t^) = 2^ exp[iq • (r - r')]Pq(r, r', uj) (102) 

where 



occ unocc 



Pq(r,r',a;) =^^ ^ Mk„(r)uk+q„'(r)uk+q„'(i-')wkn(r') 

spin kn n' 

1 1 



LU — ek+qn' + f kn + iS UJ + £k+qn' — ^kn — iS 



(103) 



The function Pq(r,r',cj) is periodic in r and r' separately and it need be calculated within a 
unit cell only which distinguishes this approach from the direct real-space method where one of 
the position variables is not restricted to the central cell. The former approach scales as N^, 
similar to localized-basis methods. It was found that the crossover between the mixed-space and 
reciprocal-space methods occurs for unit cells as small as that of Si. 

The real-space or the mixed-space approach is suitable for systems with large unit cells and a 
large variation in the electron density, or open systems. 



IV. SIMPLIFIED GW 

G W calculations are computationally expensive and therefore it is desirable to find simplifications 
which reduce the numerical effort but still maintain the accuracy of the full calculations. So far 
there is no simplified GW^ theory that is applicable to all systems. Most simplified theories are for 
semiconductors and insulators. Although some success has been achieved, none of these models 
is very reliable. While band gaps can be reasonably well reproduced by these models, details 
of bandstructures are not satisfactorily described. It is a major challenge to construct a good 
approximation for the self-energy which includes both non-locality and energy dependence but 
which is simple enough to be applicable to complex systems without significant loss of accuracy. 



A. The static Coulomb-hole and screened-exchange (COHSEX) approximation 

One of the earliest attempts to simplify G FF self-energy is the static COHSEX approximation 
(Hedin 1965a). It is obtained formally by setting E — ekn = in equations ( |70| ) and (JT^) yielding 

occ 

I]sEx(r, r') = - E 0k„(r)</)L(r')W^(r, r', c. = 0) (104) 

kn 

ScoH(r, r') = U{v - r') [W{v, r', c^ = 0) - t;(r - r')] (105) 

The first term is the exchange self-energy but with a statically screened interaction. The second 
term is the Coulomb-hole term which is the interaction energy between the quasiparticle and the 
potential due the Coulomb hole around the quasiparticle, as a result of the rearrangement of the 
electrons (screening). Both the screened-exchange and the Coulomb- hole terms in equations ( [TOI ) 
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and (|7l| ) are energy dependent and non-local respectively but in the static COHSEX approximation 
they are energy independent and the Coulomb-hole term is local. The validity of the COHSEX 
approximation relies on whether E — e^n is small compared to the energy of the main excitation 
in the screened interaction which is essentially the plasmon energy. Comparison to the results 
of the full calculations shows that the COHSEX approximation consistently overestimates the 
magnitude of the self-energy, by about 20 % in Si, resulting in larger band gaps in semiconductors 
(Hybertsen and Louie 1986). For E corresponding to an occupied state, most of the error resides 
in the Coulomb-hole term. The approximation E — e^n = is more severe for the Coulomb-hole 
term than for the screened-exchange term because the Coulomb-hole term involves a sum over 
unoccupied states as well as over occupied states whereas the screened-exchange term involves a 
sum over occupied states only. 

Another source of error comes from the neglect of dynamical renormalization or Z-factor in 
equation (p7|). When taking the matrix element (0kn|S|0kn) and setting -E — Ekn — 0, this approx- 
imately corresponds to calculating the self-energy at ii^ = Skn but the true self-energy should be 
calculated at the quasiparticle energy, leading to the formula in equation (^^. Indeed, a significant 
improvement can be obtained if the Z-factor is taken into account when the self-energy itself is 
calculated within the COHSEX approximation. 



B. Improving the COHSEX approximation 

An attempt to include dynamical renormalization into a simplified GW scheme was made by 
Bechstedt et al (1992). In this scheme, the energy dependent part of the self-energy is expanded 
to linear order around the LDA eigenvalue. A proper calculation of the self-energy derivative, or 
equivalently the Z-factor, requires unfortunately a full frequency dependent response function. It 
is not sufficient to calculate the static response function and its energy derivative at w = 0. For 
the screened-exchange term, on the other hand, knowledge of the first energy derivative of the 
response function at w = is sufficient to expand Ssex to first order around the LDA eigenvalue. 

To take into account the energy dependence of the self-energy, the screened interaction W is 
calculated within the plasmon-pole model. The static dielectric function is modelled by (Bechstedt 
et al 1992) 



e(q,p) = 1 



(.-D-. + .f^V.?/-^ 



\QTF J 4 ykpqTF 



(106) 



where kp and qxp are the Fermi and Thomas-Fermi wave vectors respectively which depend on 
the average electron density p. This model for e interpolates between the free-electron gas result 
1 -|- {2uip/q^Y at high (?, 1 -t- a{(jTF I'if' at small q (Thomas-Fermi theory) and the (7 = value 
eo for the semiconductor. It allows analytical calculation of the screened Coulomb interaction. 
The coefficient a is obtained by fitting the model dielectric function to a full RPA calculation for 
w = 0. The values of a turn out to be material independent for those semiconductors considered 
(Si, GaAs, AlAs, and ZnSe). 

To take into account local field effects due to the inhomogeneity in the charge density, an LDA 
ansatz is used (Hybertsen and Louie 1988a) 

W(v, r', w = 0) = i \W'\v - r', p(r)) + VF'*(r - r', p(r'))] (107) 

where W^ is the screened interaction for the electron gas. The sum rule, that the total induced 
charge around a test charge is —1 + 1/eo, is fulfilled but the induced charge density is allowed 
to vary according to the local density. Using this model, the static Coulomb-hole term can be 
calculated analytically (Bechstedt et al 1992) 



T-coH(r) = -(l-lV'^--('') 



eo / V" 



qTF (r) / 3eo 



Ykp (r) V eo - 1 



(108) 
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It is interesting to observe that the matrix element of S''^" = E — Scohsex and its energy 
derivative can be shown to be independent of the state if local field effects are neglected. The state 
dependence thus comes from the local field effects and it is found to be rather weak. We note 
also the dependence on qtf ~ P^ as in the LDA. E''*'" gives an upward shift of ^1.4 eV for all 
the materials mentioned above. This means that good values for the band gaps are obtained by 
using S = EcoHSEX but taking into account the dynamical renormalization factor Z. The results 
for Si, GaAs, AlAs, and ZnSe show agreement with the full Gl^F calculations to within 0.2 eV for 
most of the states considered. Application was also made to GaN recently (Palummo et al 1995). 
This scheme reduces the computational effort by two orders of magnitude. Further tests on a 
wide range of semiconductors and insulators are desirable to evaluate the validity of the number 
of approximations used in the model. 



C. Extreme tight-binding models 

One of the major problems in GFF calculations is the calculation of the response function. In 
electron-gas-like materials such as the alkalis, it is reasonable to model the dielectric function with 
a single plasmon and to neglect off-diagonal elements. In semiconductors, however, local field 
effects, which are described by the off-diagonal elements of the dielectric matrix, are important. 
The local field is dominant at distances on the interatomic scale so that high Fourier components 
are needed to describe it in plane wave basis resulting in a large matrix. It is therefore advantageous 
to calculate the dielectric matrix and its inverse directly in real space. Real-space approach is also 
useful for systems with low symmetry. Ortuno and Inkson (1978) proposed an extreme tight- 
binding model where the valence and conduction bands were assumed to be flat so that the band 
gap Eg is the only parameter entering the model. Due to its simplicity, the model allows analytic 
evaluation of the dielectric function to a certain degree. 

The wavefunctions are expanded in localized Wannier functions 

V'kn (r) - ^ ^ e^'^-'^c^nu (r - T) c„, (k) (109) 



where (j)ni> is a Wannier function localized in bond v and T is a lattice translational vector. In the 
two-flat-band model, the coefficients Cnv are equal to one. Since (pnv is localized in each cell and 
in each bond, we have 

/ d?r (bn, (r - T) 0„^ (r - T') ex 5,^5tt' (HO) 

Using this expansion in the expression for the polarization function and using the flat-band ap- 
proximation, eic'n' — eicn = Eg, one arrives at a simple expression for the dielectric function: 

e (r, r', lu) = 5 {r - r') - N (lo) ^ / d^r"v (r - r") A^ (r" - T) A* (r' - T) (111) 



where 

4K 
NH= , ,' .^„ (112) 

W^ — {Eg — 10)'^ 

A, (r - T) = 0eo„d.. (r - T) cj^l^,^^ (r - T) (113) 

This has the form of a separable matrix e = 1 — BC (Hayashi and Shimizu 1969, Sinha 1969) 
whose inverse is given by e^^ = 1 + B (1 — CB) C. Under certain approximations, the matrix e 
can be inverted giving a screened interaction (Ortuno and Inkson 1979) 
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w 



{v,v',uj) = /'dVe-i(r,r",c^)t.(r"-r') 



= w (r — r') — 



^E„ 



El+u^l-u^ .T 



Y,D,{v-T)Dl{v'-T) 



(114) 



where 



D^{v-T)= I d^r'v (r - r') A^ (r' - T) 



(115) 



and LOp is the plasmon energy. The approximation is good in the hmit Eg ^ uip. The screened 



potential is effective at uj = Eg and approaches a bare value at oj 



The quantity D^ (r - T) 



represents a dipole moment at r due to bond v in cell T. The physical interpretation of the second 
term in equation (114) is that an electron on a site interacts with other electrons through the 
Coulomb interaction, inducing dipole moments on the other sites. These dipoles in turn produce 
a potential which is screened by other induced dipoles by the frequency-dependent factor arising 
from the inversion of the dielectric function (Sterne and Inkson 1984). 

Using the above screened interaction, the self-energy can be evaluated analytically. Defining a 
state-dependent local potential 



y'dVl](r,r', 



Lo = E„) ^k„ (r') = y„^>k„ (r) 



(116) 



gives two potentials for the valence and conduction states respectively (Sterne and Inkson 1984) 

.11/3 



Vr 



COTSEX (r) = -p'^' (r) 



T/cond 
^COHSEX 



ir) = -p'/Hr)l 
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where 



C = 



(frcPr'Al (r - T) w (r - r') A^> (r' - T) 



is an on-site exchange interaction and 



eo 



El 



n 



(119) 



(120) 



Comparison with the LDA exchange-only potential, V^f^^ = — [S/vr] ' p^''^, gives a value for 

7 = 2 [9/ (27r^)] . It is interesting to observe that although the potentials have been derived 
from extreme tight-binding picture, one arrives at a formula similar to that obtained from the 
electron gas (Sterne and Inkson 1984). 

The method has been applied to C, Si, GaAs, Ge, and ZnSe (Jenkins, Srivastava, and Inkson 
1993). Good agreement with the full GW calculations are obtained for Si and GaAs but not 
for the other two materials. Perhaps this is not surprising in view of the flat extremal bands in 
Si and GaAs, making them ideal for the present approach. Moreover, the pseudopotentials used 
in the LDA calculations are better for the first two materials than for the last two so that the 
eigenf unctions used to calculate the matrix elements of the self-energy are correspondingly better. 
The LDA exchange only gap is 0.41 eV whereas the calculated value is 1.25 eV (experimentally 
1.17 eV). The direct gap at the F point is also improved from 2.49 eV to 3.32 eV (experimentally 
3.35 eV). In the case of Ge the improvement is not as good as in the case of Si but while the 



27 



LDA predicts Ge to be a metal with a negative band gap of -0.19 eV, the method gives a gap of 
0.34 eV albeit too low compared to experiment (0.89 eV). In diamond, the method improves the 
gap at the F point from 5.46 eV in the LDA (exchange only) to 7.87 eV (experimentally 7.4 eV). 
Similar improvement is also obtained across the Brillouin zone. Finally, the results for ZnSe are 
rather poor. The band gap is overestimated by ~ 1 eV (3.77 eV vs 2.82 eV). The problem could 
be related to the difhculties of treating 3d states within the pseudopotential approach resulting in 
relatively poor wavefunctions. 

A similar tight-binding approach was also proposed by Hanke and Sham (1988) and Bechstedt 
and Del Sole (1988). They derived an analytical model for E, V™, and the gap correction in 
insulators. They arrived at a similar p^ " formula for the exchange-correlation potential but they 
also included a term corresponding to the bare conduction-band exchange which is neglected in the 
Sterne and Inkson model. The method was applied to Si, diamond, and LiCl and the results are 
within a few percent of the more elaborate calculations. A more recent work uses orthogonalized 
linear combination of atomic orbitals with applications to diamond, Si, Ge, GaAs, GaP, and ZnSe. 
The gaps are generally within 10 % of the experimental values (Gu and Ching 1994). 



D. Quasiparticle local density approximation (QPLDA) 

This approximation is based on the work of Sham and Kohn (Sham and Kohn 1966) who 
showed that for a system with a slowly varying density the self-energy is short range in |r — r'| 
although for larger energies the self-energy is expected to be long range since screening is not very 
effective at large energies. Numerical calculations by Hedin (1965a) show that the self-energy of 
the homogeneous electron gas quickly vanishes beyond |r — r'| = 2rs for lu sa Ep and it depends 
therefore on the density in the vicinity of (r + r') /2 which suggests the following approximation 
(Sham and Kohn 1965) 

Eir,r',Lu)^j:hir~r',u;~Afi{n);n) (121) 

where T^h (r — r', E; n) is the self-energy operator of the homogeneous electron gas with density n, 
n = n [(r -I- r') /2] , and A/i = fx — ^h {n) is the difference between the true chemical potential and 
that of the homogeneous electron gas of density TT The local density approximation for the self- 
energy is obtained by assuming the quasiparticle wavefunction in equation (Eq) as a superposition 
of locally plane- wave-like functions (Sham and Kohn 1965) 

* (r,tj) w A{r)e'^^''^'^'>-'' (122) 

which when inserted into the quasiparticle equation gives the solution lu (E) = E ii k satisfies the 
local-density condition k = k^u with 

-E+ ^klo + V" (r) + S {kLD,E - Afi (n) ;n) = (123) 

In obtaining the above equation, the r dependence of A and k has been neglected. For a slowly 
varying density, the local change in the chemical potential from its homogeneous value is simply 
given by the electrostatic potential (Thomas- Fermi) : 

V" (r) = A/i [n (r)] (124) 

and by definition 

f^h{n) =-4(n)-f/i"^(n) 

-i4(n) + S(fcF,/iftH,«) (125) 

Using the last two expressions, the condition for the local-density wave vector becomes 
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- {kio - ki) ={E-^i)- [^h {kLD, E-^^i (n) ; n) - E^ (fc^, M/. (") ; n)] (126) 

When operating on the locaUy plane-wave function, the effect of the non-local operator 
Eft,(r — r',i5) can be reproduced exactly by a local operator Yjh{kLD,E\n)5 {v — y') . The lo- 
cal density approximation for the self-energy operator consists of using this local operator when 
operating on the actual wave function at energy E. It may happen that for some values of E 
equation ( |126| ) has no solution for positive k^j^. In this case, one should analytically continue the 
self-energy operator into complex momentum. An early work of approximating the self-energy by 
a local potential using information from the local density is given by Hedin and Lundqvist (1971). 
Instead of calculating the full self-energy E/j {k,E-,n) it is useful to write (Wang and Pickett 
1983, Pickett and Wang 1984) 

E/i (fc, E; n) = ^i"^" (n) + E^ (fc, E; n) - E^ [kp, ^J,h^, n) 

= Ai"= {n) + A (fc, E- n) (127) 

since fj^'^ has been calculated rather accurately for the electron gas. The GWA is then used to 
calculate the relatively small quantity A. Hedin and Lundqvist (1971) calculated A(fcLDA) and for 
metals, the self-energy corrections in the local-density approximation turn out to be very small. 
A is of the order of a few xlO^^ eV. For Ni, for example, the corrections are smaller than 0.1 eV 
(Watson et al 1976). The main reason for such small corrections is probably due to the weak energy 
dependence of the electron-gas self-energy. Furthermore, most of the effect of the self-energy is 
already accounted for by the local term ji'^'^. The QPLDA self-energy is dominated by a constant 
contribution and the assumption of slowly varying density is probably not good enough in many 
cases. 

Above, the homogeneous electron gas was used to obtain an approximate self-energy for inho- 
mogeneous systems. For semiconductors and insulators one can also use a homogeneous insulator 
model for this purpose which turns out to give a significant improvement. The presence of an 
energy gap in these systems brings in new physics not presence in the electron gas. The response 
function needed to calculate the screened interaction is obtained from a model semiconducting 
homogeneous electron gas. One simply introduces a gap in the otherwise metallic spectrum of the 
electron gas, that is to say (Penn 1962, Levine and Louie 1982) 

/, N f 0, \uj\ < XEp .,„c^ 

''^''^^^^{erHk,^-),H>XEp (128) 



where a;_ — y^uj^ — X^Ep sgn uj. This function satisfies the /-sum rule. The dielectric gap XEp, 
which is in general larger than the minimum direct gap, is determined from the experimental value 
of the static, long-wavelength dielectric constant eg, 

XEf = ,^^ (129) 

where Wp = v47m is the plasmon frequency. The single-particle spectrum has a gap Eg 

^(k) = ifc2 + i^gSgn(A;-fcf) (130) 

The gap Eg is taken to be the average value across the Brillouin zone of the direct gap and it 
is numerically different from XEp. It can be shown that the valence bands are lowered and the 
conduction bands are raised by the self-energy corrections, independently of the density sampled 
by the wave functions. This means, the QPLDA will always lead to an increase in the gap over 
the LDA value (Wang and Pickett 1983, Pickett and Wang 1984). 

The QPLDA was applied to Si and diamond with considerable success. For Si, the zone-boundary 
transitions A"4 ^ ATi, Lg -^ L3, and Lg — > L\ are well reproduced by the QPLDA. The zone-centre 
transition r25' — > Pis is underestimated by 0.25 eV, which is the worst case but the transition 
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r25' -^ Tis is reproduced well. The calculated indirect gap of 0.93 eV still underestimates the 
experimental value by 0.24 eV. The QPLDA gives a valence-band width of 12.5 eV, in good 
agreement with the experimental data 12.4 ± 0.6 eV (Grobman and Eastman 1972) and 12.6 ± 
0.6 eV (Ley et al. 1972). For diamond, the QPLDA gives an indirect gap of 5.74 eV, which is 
a slight overestimate compared to the experimental value of 5.47 eV. The LDA value is 4.05 eV. 
The r25' -^ ri5 gap is experimentally about 7.3 eV and the QPLDA gives 7.36 eV (compared 
to the LDA value of 5.51 eV). The valence-band width is 23.4 eV, in good agreement with the 
experimental data of McFeely et al. (24.2 ±1.0 eV) (McFeely et al. 1974) but in disagreement 
with the data of Himpsel et al (21 eV) (Himpscl, van der Veen, and Eastman 1980). The LDA 
band width is 20.4 eV. 



E. LDA + 5EcoHSEx 

In metals, the screened interaction W decreases rapidly for |r — r'| > k^\, but in insulators and 
semiconductors, it decreases as 1/eo |r — r'| for large |r — r'| since the screening is not complete. 
Accordingly, one writes (Gygi and Baldereschi 1989) 

W (r, r', Lo) ^ W^^'^ (r, r', uj) + SW (r, r', u) (131) 

where W^^*^ is the short-range interaction potential of a metallic inhomogcneous electron gas and 
SW has the same long-range behaviour as W. The self-energy arising from the first term in W has 
been shown by Sham and Kohn (1966) to be short ranged and depends only on the density in the 
vicinity of r and r', and is therefore approximated by a local potential: 



■^GW 



(r, r', iv) = fi""" (r, uj)5{r-r') + — f duj'G (r, r', u + lo') 5W (r, r', u') (132) 

2tt J 

Pickett and Wang (1984) found that the inclusion of energy dependence in the local exchange- 
correlation potential had very little effects on the eigenvalues obtained with an energy-independent 
LDA potential. It is then assumed that for states close to the Fermi level, jx^'^ takes its value at 
the Fermi level nYhh ('') • Furthermore, it is assumed that 5W depends only on |r — r'| , which is 
strictly valid only when |r — r'| -^ oo, but it has been found numerically that in non-metals the 
local field effects become negligible as |r — r'| exceeds the interatomic distance (Hybertsen and 
Louie 1986, Godby, Schhiter, and Sham 1988). One now makes the COHSEX approximation on 
i5S, rather than on the full self-energy: 

5^ (r, r') - -p (r, r') 5W (r - r') + ,5Scoh (133) 

where the first term is the screened exchange contribution and 

5ScoH - — ^ / d\5W (q) (134) 

2 (27r) J 

is the Coulomb hole contribution which is a constant, shifting all eigenvalues by the same amount. 
The Fourier transform of 5W is given by 

5W (q) = ^ [e-'c (q^ q^^ - 0) - 6^1 (q, u; = 0)] (135) 

where e^^ (q, q,0) is the diagonal part of the inverse dielectric matrix for the semiconductor calcu- 
lated in the RPA and ej]^{q,0) is the inverse of the static Lindhard dielectric function (Lindhard 
1954) of a homogeneous electron gas. 

Applications to diamond. Si, Ge, GaAs, AlAs, and GaP give in general agreement to within 0.2 
eV with the results of the full calculations. Larger deviations (0.3-0.4 eV) occur for the X4y -^ Xic 
gap in diamond and the T25'v ^ ^2'c gaps in Ge. The Coulomb-hole contribution SEcoR is positive 
and cancels approximately the effect of the screened exchange term in the valence bands. As a 
result, the net effect of the self-energy correction is to shift the conduction band upwards and leave 
the valence bands essentially unchanged. This is also in agreement with the full calculations. 
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V. APPLICATIONS 

A. Alkali metals 

The early calculations in the GWA were performed for the electron gas, because of its simplicity. 
(For a review see Hedin and Lundqvist 1969). The alkali metals are therefore of particular interest, 
being the systems which are the closest approximation to the electron gas. For these systems the 
correlation effects are only moderately strong, and the GWA could therefore be expected to be 
relatively accurate. It therefore created a lot of interest when the first high-resolution angular 
resolved photoemission experiments for the alkali metals appeared (Jensen and Plummer 1985, 
Lyo and Plummer 1988, Itchkawitz et al 1990). Two unexpected observations were made. First 
the band width was smaller than expected. For Na the occupied part of the band was found to 
be 2.65 eV broad (Lyo and Plummer 1988). This is consistent with an old angular-integrated 
measurement by Kowalczyk et al (1973). This band width is substantially smaller than the width 
3.2 eV predicted by nearly free electron (NFE) theory. Figure shows the experimental peak 
position (crosses) as a function of the photon energy together with the prediction of NFE theory 
(dashed line). The band narrowing is about a factor of two larger than the narrowing (0.27 
eV) predicted by GW calculations for the electron gas of the appropriate density (Hedin 1965a). 
Although this error is not very large in absolute terms, it raised questions about the accuracy of 
the GWA for these systems. In addition a set of essentially dispersionless states were observed 
close to the Fermi energy, which was completely unexpected according to the NFE band structure. 
This raised the interesting issue that the alkali metals may have charge density waves (Overhauser 
1985). 

A GW calculation was performed by Northrup et al (1987, 1989) for Na metal and by Surh et al 
(1988) for K, using a generalized plasmon pole approximation. These calculations gave essentially 
the same band narrowing (0.31 eV) for Na as a GW calculation for the electron gas. A much 
better agreement with experiment was, however, obtained by using an improved dielectric function 
(Northrup et al 1987, 1989) 

e-i = l + i;[l-P(t; + X^^)]-ip, (136) 

where P is the independent particle polarizability and X™ — 5V^'^/5p with V^'^ the (LDA) 



exchange-correlation potential and p the density. Equation (136) gives the appropriate dielec- 
tric function within the density functional formalism. Using this e instead of the RPA e, the band 
width narrowing increased from 0.31 cV to 0.57 eV. Finally, a calculation was performed using a 
Green function with a certain self-consistency, namely with the LDA eigenvalues replaced by the 
quasi-particle energies. This led to a further increase of band width narrowing to 0.71 eV. Lyo and 
Plummer (1 988) also observed large effects of including the corrections to the dielectric function 
in equation ( |136| ). 

The theoretical justification for including such corrections without simultaneously adding vertex 
corrections is, however, weak. We notice that the total energy of the system can be expressed in 
terms of the dielectric function. The quasi-particle energies can then be obtained by differentiating 
with resp ect t o the occupation numbers (Rice 1965). If the dielectric function has the form of 
equation ( |136| ), it was shown that there is also a vertex correction of the type (Rice 1965, Ting 
1975, Mahan 1994) 

r = T^^. (137) 

Mahan and Sernelius (1989) have extensively tested various corrections to the dielectric function, 
including the corresponding vertex corrections when calculating the self-energy. They found that 
for the electron gas the vertex corrections cancel most of the effects of the corrections of the 
dielectric function, and the final results are rather close to the original GWA. Del Sole, Reining, 
and Godby (1994) obtained similar conclusions for Si. The results for the electron gas suggest that 
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there are discrepancies between appropriate self-energy calculations and the peak positions in the 
photoemission experiments. 
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FIG. 4. Peak position as a function of photon energy for photoemission at normal angle from the Na 
(110) surface. The crosses show experimental results (Jansen and Plummer 1985), the filled circles show 
the photoemission calculations of Shung and Mahan (1987), the full line the quasi-particle energies and 
the dashed lines the NFE theory. (After Shung and Mahan 1987). 

The relation between the peak position in photoemission and the quasiparticle energies was 
studied by Mahan and coworkers in a series of papers, focussing on surface and quasiparticle 
life-time effects (Shung and Mahan 1986, 1988, Shung ei al 1987). They studied a model which 
includes the rapidly varying potential in the surface region as well as lattice potential in the 
bulk. The effects of these potentials on the photoemission process were included, while the lattice 
part of the potential was neglected when calculating the states. The self-energy was calculated 
using the Rayleigh-Schrodinger perturbation theory, which gives results that differ slightly from 
the traditional GW results. For instance, the Na band width is reduced by 0.37 eV due to this 
self-energy. Mahan and coworkers furthermore included the effects of the imaginary part of the 
self-energy of the emitted electron, by using wave function (^> for these electrons which decayed 
exponentially inside the surface. Finally, the self-energy of the initial state was included as well. 
Due to the exponential decay of (^^ , the electron momentum perpendicular to the surface is not 
conserved, and there are nonvertical (non fcj^-conserving) transitions. 

Mahan and coworkers found that the broadening of the initial state and the instrumental reso- 
lution as well as the interference between bulk and surface photoemission shift the apparent peak 
positions towards lower binding energies relative to the quasiparticle energies by about 0.2-0.4 eV 
(Shung and Mahan 1988). Including these effects as well as the quasiparticle self-energy shifts leads 
to the filled circles in figure 0. There is a substantial band narrowing relative to the quasiparticle 
energies (full line) and the agreement with experiment is generally quite good. This suggests that 
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the GWA gives good quasiparticle energies for Na, and that the main reason for the discrepancy 
between these energies and the photoemission peak positions can be explained by considering the 
details of the photoemission process. 

Mahan and coworkers also gave an explanation for the unexpected features close to the Fermi 
energy (see the essentially dispersionless structure in figure Q). These features occur for photon 
energies for which there are now vertical, energy-conserving transitions available. Due to the 
exponential decay of (fr' , nonvertical but energy-conserving transitions are available. It was shown 
that interference between surface and bulk emission puts most of the weight of these transitions 
close to the Fermi energy (Shung and Mahan 1986). figure || illustrates that the theory can almost 
completely explain the experimental features, with just a few experimental points at larger photon 
energy unexplained. Thus there seems to be no need to assume a charge-density wave to explain 
this structure. 



B. Semiconductors and Insulators: sp systems 

The LDA systematically underestimates the band gaps in semiconductors and insulators. In 
table I the calculated LDA band gaps of some materials are compared with the experimental gaps. 
The discrepancies range from 30-100 % and for Ge the LDA conduction and valence bands in fact 
overlap when relativistic corrections are included. Also, individual bands away from the Fermi 
level can be in error up to 50 %. It is natural to ask if the band-gap problem originates from the 
error in the LDA. Exchange-correlation potentials V^^ calculated from GW self-energies turn out 
to be similar to the LDA y™ (Godby, Schliiter, and Sham 1988) which indicates that even the 
exact V'^'^ probably does not give the correct gap but this is still an open question. The V^'^ may 
be a non-analytic function of the particle number (Almbladh and von Barth 1985b, Perdew and 
Levy 1983, Sham and Schliiter 1983). That is to say, y™ with an extra electron, V^+i, may have 
an additional constant compared to V^'^ and this constant may be large. Moreover, experience 
with empirical potentials shows that a local potential cannot in general give both the correct band 
structure and the ground-state electron density (Kane 1971). 

In table I the band gaps of some materials calculated within the GWA are shown. The agreement 
with experiment is very good, in most cases to within 0.1 eV. Although it is not strictly true, the 
self-energy correction is approximately an upward rigid shift of the conduction band relative to the 
valence band, the so called scissor operator, i.e. cutting the band structure along the band gap 
and shifting the conduction band rigidly upwards. The scissor operator is accurate to 0.1, 0.2, 0.2, 
and 0.4 eV in Si, GaAs, AlAs, and diamond respectively (Godby, Schliiter, and Sham 1988). The 
validity of the scissor operator in Si is somewhat fortuitous, due to an almost complete cancellation 
between the strong energy dependence and non-locality. 

The experimental values are obtained from optical measurements or photoemission and inverse 
photoemission experiments. The latter corresponds closer to the theoretical values whereas the 
former may contain excitonic binding energy which should be subtracted off but unfortunately it is 
unknown. In optical experiments, the excited electron does not leave the system and may therefore 
form an exciton with the corresponding hole. 

To discuss in more details the features in the self-energy which are important for the quasi- 
particle energies, we consider Si as a prototype since it has been studied extensively. The main 
features are energy dependence and non-locality. We first consider non-locality within the COH- 
SEX approximation. A measure of non- locality in the self-energy is its range, the distance |r — r'| 
beyond which the self-energy is approximately zero. This range r^ is approximately given by the 
corresponding value for a jellium with the average density of Si (r^ = 2) and r^ ^ 'irg. In Si 
more than 99 % of the matrix element of the self-energy in a state ip, {ip\Y,\ip), originates from 
|r — r'l < r/j. One expects non-locality to be important when rh is comparable to or greater than 
the extent or wavelength of the wavefunction which is the case in Si. The matrix element in a 
non-local potential can be very sensitive to the nodal structure of ip. This is contrary to the case 
when E is local such as V^'^ since it is then |-0p that enters into the integral. In Si, non-locality 
has the effect of widening the gap. This can be understood as follows. The top of the valence band 
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is bonding p whereas the bottom of the conduction band is antibonding p. Therefore non-locality 
has a larger effect on the conduction band than on the valence band since the antibonding state 
has an extra node which means a smaller wavelength than that of the bonding state. The presence 
of an extra node in the antibonding state reduces the matrix element (V'|S|')/') relative to {ilj\V^'^\ip) 
and therefore the conduction band is pushed upwards. Non-locality has a smaller effect on the 
valence state so the net effect is a widening of the gap (Godby, Schliiter, and Sham 1988). 
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FIG. 5. (a) and (c) show the screening potential in response to a single electron at r' (indicated by -f ) 
in the (110) plane of Si in units of Ry. (b) and (d) show the contribution from local fields only. The bond 
chain is indicated by a straight line. After Hybertsen and Louie (1986). 



The non-locality arises from the density matrix (exchange charge) 

occ 

p,(r,r') = ^V.(r)V:(r') 



(138) 



and the screened potential T4^(r,r',u;). In semiconductors, the screening is not complete because 
due to the gap, a finite energy is required to excite a particle-hole pair. The range of the screened 
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interaction and its nodal structure are then determined by Px ■ Apart from non-locality, anisotropy 
also plays an important role. In a homogeneous system, px and W only depend on |r — r'|. In an 
inhomogeneous system, the screening potential W — v a.s a. function of r' may strongly depend on 
the location r of the test charge. In Si, for example, there is a large accumulation of charge in 
the bonding region as opposed to the antibonding region. It is to be expected that the screening 
potential of a test charge located in the bonding and antibonding regions will be very different, 
as shown in figure ^ This local field effect , which is entirely missing in the homogeneous case, is 
very important in covalent materials. Local fields are crucial in determining the strength of the 
screening hole but not its shape and they contribute directly to the differing strengths of S at 
different points in the unit cell and therefore to the band-gap correction figure H (Hybertsen and 
Louie 1986). In Si local fields account for more than one-third of the screening potential in the 
region around the bond as can be seen in figure ^. 




FIG. 6. Contour plots of the self-energy E(r, r', (j;=inidgap) in eV a.u.~'^ for r fixed at the bond centre 
and r' shown in the (110) plane for (a) Si, (b)GaAs, (c) AlAs, and (d) diamond. For Si the corresponding 
plots with r fixed at the tetrahedral interstitial site is also shown (e). For comparison, the self-energy 
operator of jellium with Ts — 2.0 (the average electron density of Si) is shown in (f) (from Hedin (1965a)). 
After Godby, Schliiter, and Sham (1988). 

The local fields shift the centre of the screening potential and increase it in the bonding region 
and reduce it in the interstitial. The local field effect is one important feature in Si which distin- 
guishes covalent-bond semiconductors from the alkalis. The screening potential in these materials 
is therefore considerably more complex than in the alkalis. Calculations for Si show that the local 
field effect in W is confined to a region near r as shown in figure o. This is not surprising since 
the long-range part of the potential is contained in the diagonal element VF(q) corresponding to 
small q. In the plane-wave basis, the local fields are described by the off-diagonal elements of the 
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dielectric matrix. Nevertheless, Y,{r,r',uj = 0) is almost spherical as a function of r' for a fixed 
r and it is reasonably well reproduced by jellium S with the same average density as that of Si. 
This is illustrated in figure |6[ However, the interaction of the wavefunction with the non-locality 
in S is not contained within the jellium model. 

The local field effect has a large influence on Ecoh- In fact, in a homogeneous system Scoh 
is a constant within the COHSEX approximation (equation ( |105| )). Thus, if local field effect is 
neglected, Scoh has no effects on the band dispersion. Scoh is deeper in the bonding region and 
shallower in the antibonding region as shown in figure M. Since the valence state is concentrated 
in the bonding region and the conduction state in the antibonding region, Scoh makes a large 
contribution to the gap. Scoh is, however, a local potential within the COHSEX approximation, 
a feature which is common to V^'^ as opposed to Ssex which is non-local. This mechanism of gap 
opening by Scoh could in principle be accounted for by V^"^. The local field effect on Ssex, on 
the other hand, is small, since the local field effect is rather localized and Ssex is dominated by 
the large bare Coulomb interaction for small |r — r'|. The local field contribution to Ssex is about 
25 % of that to EcoH but of opposite sign (Hybertsen and Louie f986). Non-locality is essential 
in determining the correct quasiparticle energies and in particular the band gap. 




V COH ^ -' 



FIG. 7. The local potential corresponding to the Coulomb hole part of the electron self-energy in the 
COHSEX approximation in the (110) plane of Si in units of Ry. After Hybertsen and Louie (1986). 

The COHSEX approximation is a static approximation to the self-energy. Strictly speaking, the 
quasiparticle energy should be obtained fronr the position of the peak in the spectral function. 
This procedure requires knowledge of the energy dependence of the self-energy at least around the 
quasiparticle energy. The degree of the energy dependence is related to the renormalization factor 
Z (weight of the quasiparticle, equation (p7|)): the smaller Z, the stronger the energy dependence. 
The value of Z for semiconductors is 0.7-0.8 (Hybertsen and Louie 1986). The slope of the real 
part of the self-energy around the chemical potential is negative and the renormalization factor Z 
reduces the self-energy correction, as shown in equation (^ . Thus, if the energy dependence of the 
self-energy is neglected, that is if the self-energy were calculated at the LDA instead of quasiparticle 
energy, the band-gap correction would be overestimated. This is in agreement with the results of 
the COHSEX approximation which approximately corresponds to neglecting the renormalization 
factor Z. A similar conclusion is reached if the self-energy is approximated by its value at w = 0. 
In this case the self-energy correction for the valence state would be underestimated whereas for 
the conduction state overestimated, leading again to an overestimated band gap. A much better 
agreement is obtained if the self-energy correction (Ecohsex — ^™) is simply multiplied by Z. We 
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note that the COHSEX approximation without local field effect generally gives a gap in better 
agreement with experiment, although not for Si. This means that local field effect, or in a certain 
sense non-locality, tends to cancel energy dependence. The importance of energy dependence is 
illustrated by plotting (V'kn|S(£'kn) — V'™|'(/'kn) and (V'kn|S(0) — V^'^'^IV'kn) as a function of Skn- 
The first quantity is very much like a step function (scissor operator) while the second quantity 
shows a strong energy dependence (Figure 5b and 6b of Godby, Schliiter, and Sham 1988). The 
effect of energy dependence is therefore to alter greatly the dispersion of the individual bands. 

In general, the energy dependence of S leads to a strongly state-dependent self-energy correction 
AE = E — V^'^ within each band as well as across the gap. The weak state dependence of AE 
for Si within a band, resulting in a scissor operator, is therefore a coincidence. In diamond, for 
instance, the scissor operator approximation is not as good as in Si. 
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FIG. 8. (a) and (c): Tlie excliange charge in tlie (110) plane of Si with r' fixed (indicated by -I-) in units 
of electrons/cell, (b) and (d): The screened- exchange part of the self-energy x|r — r'| in the COHSEX 
approximation in units of a.u. Ry./cell and the contours increase in powers of 2. After Hybertsen and 
Louie (1986). 

The highest occupied state in the exact DFT gives the exact ionization energy (Almbladh and 
von Barth 1985a). Assuming that the LDA y^ is close to the exact one, we expect the self-energy 
correction to shift the conduction band but not the top of the valence band. If the RPA is used 
to obtain the electron gas data (von Barth- Hcdin potential 1972), it is indeed found that the the 
top of the valence band is almost the same within the LDA and the GWA (Godby, Schliiter, and 
Sham 1988). More accurate electron gas data obtained from GW calculations (Lundqvist and 
Samathiyakanit 1969) or Quantum Monte Carlo simulations (Ceperley and Alder 1980), however, 
shift the top of the valence band upward by about 0,5 eV. Thus, vertex corrections (corrections 
beyond the GWA) may shift the GW bandstructure by 0.5 eV upwards (Godby, Schliiter, and 
Sham 1988) or the LDA may not be sufficiently accurate. 

Other authors have also repeated GW calculations for Si band structure with similar results 
(Hamada, Hwang, and Freeman 1990, Rohlfing, Kriiger, and PoUman 1993). Calculations of the 
quasiparticle bandstructure and the band gap of many semiconductors and insulators have been 
performed by a number of authors. Quasiparticle band structures of six II- VI compounds ZnS, 
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ZnSc, ZnTe, CdS, CdSe, and CdTe have been calculated by Zakharov et al (1994). We list in 
table I the results of these calculations. Quasiparticle energies for the F-center defect in LiCl have 
also been calculated by Surh, Chacham, and Louie (1995) and self-energy calculations of carrier- 
induced band-gap narrowing in Si may be found in the work of Oschlies, Godby, and Needs (1992, 
1995). Berggren and Sernelius (1981) also studied the band-gap narrowing in doped Si and Ge 
as a function of impurity concentration. Application of the GWA to metal-insulator transition 
of Si in the diamond structure (Godby and Needs 1989) suggests that the metalization occurs 
at a much smaller volume than in the LDA which in turn indicates that the the Fermi surface 
obtained from the Kohn-Sham DFT is not necessarily the same as that of the real system as 
shown by Schonhammcr and Gunnarsson (1988) for model systems. Recently, application to a 
two-dimensional crystal found good agreement between the quasiparticle energies in the GWA and 
Quantum Monte Carlo results (Engel, Kwon, and Martin 1995). 



LDA 



GWA 



Expt. 



AlAs 


1.37 


Alo.5Gao.5As 


1.12 


AIN (wurtzite) 


3.9 


AIN (zinc-blende) 


3.2 


AlP 


1.52 


AlSb 


0.99 


CdS (zinc-blende) 


1.37' 


CdS (wurtzite) 


1.36 


CdSe (zinc-blende) 


0.76 


CdSe (wurtzite) 


0.75 


CdTe (zinc-blende) 


0.80 


CdTe (wurtzite) 


0.85 


Diamond 


3.90 


GaAs 


0.67 


GaN (wurtzite) 


2.3 


GaN (zinc-blende) 


2.1 


GaP 


1.82 


GaSb 


-0.10 


Ge 


<0 


InAs 


-0.39 


Ino.53Gao.47As 


0.02 


InP 


0.57 


InSb 


-0.51 


LiCl 


6.0 


LiaO 


5.3 


MgO 


5.0 


Si 


0.52 


SiC(/3) 


1.31 


ZnS (zinc-blende) 


2.37 


ZnS (wurtzite) 


2.45 


ZnSe (zinc-blende) 


1.45 


ZnSe (wurtzite) 


1.43 


ZnTe (zinc-blende) 


1.33 


ZnTe (wurtzite) 


1.48 



O.SS'' 



2.18" 

2.06"= 

5.8"* 

4.9^* 

2.59"= 

1.64*= 

2.83\ 

2.79' 

2.01' 

1.91' 

1.76' 

1.80' 

5.6", 5.33\ 5.67*= 



2.45'' 



1.58" 

3.5** 

3.1** 

2.55" 

0.62*= 

0.75" 

0.40*= 

0.80*= 

1.44*= 

0.18*= 

9.1" 

7.4'' 

7.7-'' 

1.29" 

2.34*= 

3.98' 

4.03' 

2.84' 

2.75^ 

2.57 

2.67 



1.32^ 1.22*= 



0.65*= 



1.24", 1.25*= 



2.32*" 
2.09*" 
6.2" 

2.50*" 

1.68*" 

2.55"'" 

2.59*" 

1.90*" 

1.97*" 

1.92*" 

1.60*" 

5.48*" 

1.52** 

3.5*" 

3.2^ 3.3' 

2.39*" 

0.80*" 

0.744*" 

0.41*" 

0.8T" 

1.42** 

0.23*" 

9.4*" 

6.6* 

7.83*^ 

1.17'' 

2.39** 

3.80*" 

3.92*" 

2.96*" 

2.87" 

2.71" 



1.63^ 



TABLE I. Minimum band gaps of semiconductors and insulators which have been calculated within 
the GWA. The energy is in eV. 
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1. Core polarization 

A disadvantage of using pseudopotentials is a difficulty of including the core states in the calcu- 
lation of the exchange potential as well as the polarization. Within a pseudopotential scheme, it 
is inevitable that core electrons are treated in an approximate manner. The core electron charge 
densities are frozen at their atomic values in the reference configuration used for constructing the 
pseudopotential. The core electrons and their potential are then eliminated so that there is no 
possibility for them to relax. Core relaxation gives rise to crystal field distortion and strong mixing 
between 3d and 3p states. These single-particle effects are small and can be included a posteriori 
by comparison with all-electron calculations. It is, however, important to take into account many- 
body effects arising from core relaxation since they can be large. In atoms with easily polarizable 
core the inclusion of core relaxation leads to an increase in the ionization energies, a contraction 
of the valence shell, a reduction of polarizabilities and oscillator strengths of the valence electrons. 

The self-energy may be broken up into three terms (Hedin 1965b, Hedin and Lundqvist 1969): 

Y, = iGcW + iGvWv+iGyvRcV (139) 

where W — Wy +vRcV. The first term is the core- valence and core-core (screened) exchange which 
is essentially the same as the bare exchange since screening is ineffective for small distance. The 
second term is the self-energy of the valence electron and the last term is the screened polarization 
potential due to the core electrons acting on the valence electrons. In pseudopotential calculations, 
the first and last term together are approximated by the LDA. The absolute contribution from 
these terms is estimated to be ~1 eV for atomic Na (Hedin and Lundqvist 1969) and solid Al 
(Arbman and von Barth 1975) but the difference S — V^"^, which is a more relevant quantity, is 
much smaller. For s-p semiconductors, this difference can be significantly larger. As a result, the 
calculated direct gaps and also the orderings and splittings of the conduction bands in Ge and 
GaAs, materials of technological interest, are in disagreement with experiment. It is important 
to get the orderings of the conduction band right since they affect the transport properties and 
life-time of thermally excited carriers. Atomic calculations in Ge estimate the error for the 4s and 
4p states to be 0.3 and 0.04 eV respectively (Hybertsen and Louie 1986). If this error is taken 
into account, the result for the band gap becomes even better. In transition metals it is crucial to 
include the core electrons in the calculations of the bare exchange since the error can be as large 
as ^1 eV. 

While the core-valence exchange can be treated straightforwardly in the Hartrec-Fock theory, 
core-valence correlation is more complicated. A good approximation for taking into account core 
polarization is provided by the core polarization potential (GPP) method used in quantum chem- 
istry (Miiller, Flesch, and Meyer 1984). The physical idea behind this method is that core polariza- 
tion functions are characterized by sharp high-frequency excitations and rather insensitive to the 
valence environment. This allows evaluation of the core polarization function for isolated atoms 
and neglect of frequency dependence. Gonsider a valence electron in the presence of a core. The 
electron polarizes the core resulting in polarization p —acV (1/r) where r is the position of the 
electron with respect to the core. From classical theory of electrostatic, the electric field arising 
from this polarization at r' is given by 

3(p-r')r' -r'2p 
E(rO= ^^ l„ (140) 

The potential experienced by another electron at r' due to core polarization is then 

T4,e(r,r') = -a,f^ (141) 

[rr ) 

assuming that both r and r' are large. For a set of cores and many valence electrons, the core 
polarization potential is 
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^^CPP = -i^aefc-fc (142) 

C 

where fc is the electric field acting on core c due to the valence charges at i and all other cores. 

fc = ES^(--)^ElF;^- (143) 

where C is a cut-ofF function, Zc' is the net charge of core c' and Rcc' = Re' — R-c- Inserting fc 



into equation (142) yields (MiiUer, Flesch, and Meyer 1984) 



^cpp = -^E-JEi^'(--)+E 



— 3— 3-C(rcOC'(rcj) 

C2 CJ 



J c'#c '"-"-cc' c',c"^c ^^cc'^^cc" J 

The first term is the CPP in atoms with a single valence electron. Without the cut-off function, 
it would diverge at small r. The cut-ofF function C is semi-empirical with a parameter related to 
the size of the core (Biermann 1943, Biermann and Liibcck 1948, Bates 1947). The first term is a 
one-particle operator which is added to the pseudopotential. The second term is the two-electron 
interaction which was qualitatively discussed above. The third term is an indirect interaction of 
a valence electron with another core. This term is repulsive and cancels the attractive potentials 
from the first and last term. The latter is a core-core interaction and, together with the third term, 
essential to make sure vanishing long-range interaction with a neutral atom. 

The effective interaction among the valence electrons is given by (Hedin and Lundqvist 1969) 

Wc = v + v'^RcV + v'^RcvJ2Rc'V + ... (145) 

c c c'^c 

Re is the full or self-consistent response function of core c given by 

vRcV = Ve-e (146) 

Thus, the screened interaction becomes 

W = e-\ 

= [1 - VKcP°] ~' Wc (147) 

where P^ is the valence RPA polarization function. Using this CPP formalism, various transition 
energies for Si, Ge, AlAs, and GaAs were calculated by Shirley, Zhu, and Louie (1992). The results 
for the fundamental bandgaps of Si, Ge, AlAs, and GaAs shown in table are in systematically 
better agreement with experiment compared with previous calculations where core relaxation ef- 
fects were taken into account within the LDA only. The correction can be as large as 0.4 eV for 
the fundamental gap in GaAs. Notable also (not shown in the table) is the correct L — T — X 
ordering of conduction-band states in Ge obtained in the CPP approach. X — L and X^c — ^ic 
splittings in GaAs are also improved (Shirley, Zhu, and Louie 1992). 

The CPP formalism is relatively easy to implement in many-body valence calculations with- 
out increasing computational cost significantly. Apart from its use within the pseudopotential 
approach, it can also be used in all-electron calculations with frozen core. 
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LDA CPP Expt. 

Si 

r8„ -^ 0.85X5C 1.29 1.16 1.17 

Ge 

Tsv -^ Ttc 0.53 0.85 0.89 

Tav -> Lee 0.58 0.73 0.744 

AlAs 

r8„ -^ Xec 2.09 2.01 2.24 

GaAs 

Tsv -^ Xec 1.02 1.42 1.52 



TABLE II. Fundamental bandgaps of Si, Ge, AlAs, and GaAs calculated within the LDA and the 
GWA including core polarization within the CPP formalism (Shirley, Zhu, and Louie 1992) compared with 
experiment (Madelung O 1984). Energies are in eV. 
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C. Transition metals 

The success of the GWA in semiconductors has encouraged apphcations to more complicated 
systems of transition metals and their compounds. From numerical point of view, transition metals 
require a new approach for calculating the response function and the self-energy. The conventional 
plane-wave basis is not suitable in this case because the localized nature of the 3d states results 
in a prohibitively large number of plane waves. From theoretical point of view, it has been argued 
that the localized nature of the 3d states makes atomic approach as more suitable for explaining 
the characteristic properties of transition metals. The presence of a Fermi surface as shown by de 
Haas-van Alphen measurements, on the other hand, strongly suggests that the itinerant character 
of the 3d electrons should be taken into account. Moreover, the band widths in transition metals 
are not too small and the ratio between the Hubbard U and the band width is of order one. This is 
crucial when we consider screening of a photoemission hole where 3d electrons from neighbouring 
cells can take part in the screening whereas such possibility is absent in the atomic case. It 
seems then that RPA-type of approach such as the GWA is meaningful for these systems. GW 
calculations for transition metals have not been extensively performed. We concentrate therefore on 
two materials Ni (Aryasetiawan 1992a, Aryasetiawan and von Barth 1992b) and NiO (Aryasetiawan 
and Gunnarsson 1995) for which full GM^ calculations have been done in some details and on MnO 
for which a model GFF calculation has been performed (Massidda et al 1995a). 



1. Nickel 

Among the transition elements, Ni is the most anomalous case in many respects. The LDA 
bandstructure deviates significantly from angle-resolved photoemission data. The occupied 3d-band 
width is 30 % smaller than that of the LDA (3.3 eV vs 4.5 eV) (Hiifner et al 1972, Himpsel, Knapp, 
and Eastman 1979) and the experimental exchange splitting is half the LDA value (0.25-0.30 eV 
vs 0.6 eV) (Eberhardt and Plummer 1980). Furthermore there is the famous 6 eV satellite below 
the Fermi level (Hiifner et al 1972, Hiifner and Wertheim 1973, Kemeny and Shevchik 1975) which 
is entirely missing in the LDA or in any single-particle theory. These discrepancies are related 
to excited-state properties. An indication that single-particle theories would have difficulties in 
describing quasiparticle properties in Ni is the fact that the two lowest atomic configurations 3d^4s 
and 3(i^4s^ are almost degenerate, differing by only 0.025 eV (Moore 1958). Another indication 
of many-body effects in Ni is the unusually large quasiparticle widths (Eberhardt and Plummer 
1980) - up to 2 eV at the bottom of the 3d band - which implies strong interaction between the 
quasiparticles and the rest of the system resulting in short life-times. The photoemission process 
introduces an additional 3d hole to an already existing one, causing on-site many-body correlations 
not amenable to a single-particle treatment. This is in contrast to Cu where there is only one 3d 
hole after photoemission and where the LDA bandstructure is good, apart from a somewhat too 
high position (0.4 eV) of the 3d band relative to the 4s band (see e.g. Jones and Gunnarsson 1989). 
Ground-state properties such as equilibrium lattice constant, bulk modulus, and magnetic moment 
are well reproduced by the LDA with the exception of the cohesive energy where the LDA value 
is about 1 eV too small (Moruzzi, Janak, and Williams 1978). 

The GFF results for Ni may be summarized as follows (Aryasetiawan 1992a, Aryasetiawan and 
von Barth 1992b): The LDA bandstructure is much improved, in particular the 3d-band width is 
narrowed by almost 1 eV, as shown in figure 0. The quasiparticle life-times are also given rather 
well by the GWA but the exchange splittings remain essentially unchanged from their LDA values 
and the 6 eV satellite is not reproduced. 

The self-energies show a number of interesting features. As an illustration, the self-energy as a 
function of frequency for the r25 state is shown in figure 11^. The imaginary part of the self-energy 
is significantly more complicated compared with those of the alkalis or semiconductors. The latter 
are typically characterized by a large peak associated with a plasmon excitation but they otherwise 
show no other distinct structures. The frequency structure of the imaginary part of the self-energy 
is determined essentially by the imaginary part of the screened interaction W ^ as may be seen in 
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equation (|6^). In alkalis or s — p semiconductors Im W is dominated by a plasmon peak which is 
then mirrored in Im E"^. Similarly in Ni, there is a strong similarity between Im W (figure El) and 
Im S'^. In transition metals, there is no well-defined plasmon excitation, rather it merges with the 
single-particle excitations forming a broad spectrum. There is a two-peak structure at about 20 
and 30 eV which is probably due to plasmon excitation. An estimate based on the electron gas 
formula gives a plasmon energy of 30.8 eV when the 3d electrons are included in the density. This 
coincides rather well with the second large peak in the two-peak structure. A smaller structure at 
around 5-6 eV originates from transitions from the occupied valence band to states just above the 
Fermi level which constitute a large density of states. Interesting to observe is the behaviour of 
Im S'^ at large frequencies. The hole and particle parts show similar behaviour and they therefore 
tend to cancel each other when one performs a Hilbert transform to obtain the real part of the 
self-energy. This justifies the use of energy cut-off in the calculation of the response function. It 
also agrees with our physical intuition that the main contribution to the self-energy should come 
from energies up to the plasmon energy. For states lying a few eV below the Fermi level, the hole 
part (negative energy) of Im S'^ has larger weight than the particle part (positive energy). This 
simply reflects the fact that the hole (occupied) states have larger overlap and correlation with 
other occupied states resulting in larger correlated part of the self-energy for the hole part. As 
we go towards the Fermi level, the hole and particle parts become of almost equal weight. This is 
the reason why the self-energy correction for states at the bottom of the band is larger than for 
those at the top of the band, resulting in band narrowing. Im E'^ around the Fermi level shows a 
quadratic Fermi liquid behaviour but it becomes linear rather quickly. 

The real part of the self-energy is obtained by Hilbert transforming the imaginary part. A 
notable feature is a large derivative at around the Fermi level which implies large renormalization 
of the quasiparticle weights. Typical values for the renormalization factor is 0.5 for the 3d states 
(Aryasetiawan 1992a). This is significantly smaller than in the electron gas (0.7) (Hedin and 
Lundqvist 1969) or semiconductors (0.8) (Hybertsen and Louie 1986) which reflects a larger loss of 
single-particle character of the quasiparticles. The s states on the other hand have renormalization 
factor '^ 0.7, comparable to those in the alkalis and semiconductors. It is in agreement with our 
physical picture that the 3d states are more correlated than the 4s-4p states. 

In comparison with semiconductors, the self-energy corrections in Ni is considerably more com- 
plicated. While in semiconductors the self-energy correction is approximately a scissor operator 
which increases the band gap by shifting the conduction band upwards and leaving the valence band 
unchanged, the self-energy correction in Ni shows a rather strong state and energy dependence. 
The self-energy correction can be positive or negative depending on the state and its magnitude 
varies throughout the Brillouin zone. For example, at the X-point the bottom of the 3d band 
experiences a self-energy correction of 0.8 eV while the top of the band is almost unchanged. The 
correction to the state T'2^ is positive whereas at L'2 state it is negative. The state dependence of 
the self-energy correction is demonstrated clearly in figure |9| by the lowest valence band which is of 
4s character at the bottom and a mixture between s and d at the top. The self-energy correction is 
approximately zero at the bottom of the band and positive at the band edges. The free-electron- 
like s states are well described by the LDA but the description of the more localized d states is 
less satisfactory. 

As can be seen in figure O, the 3d band width is reduced by almost 1 eV. The exchange splittings 
on the other hand remain essentially unchanged from their LDA values. The discrepancy between 
the LDA exchange splitting and the experiment is rather small, 0.3 eV, which is slightly larger 
than the numerical accuracy (0.1-0.2 eV) but the results seem to indicate inadequacy in the GWA 
itself. G FK calculations on transition metal atoms also show that strong correlations among 3d 
electrons of opposite spin are not well accounted for by the GWA (Shirley and Martin 1993). The 
quasiparticle widths or the inverse life-times are in reasonable agreement with available experimen- 
tal data (Aryasetiawan 1992a). The large width at the bottom of the 3d band indicates a strong 
interaction between the quasiparticles and the rest of the system. 
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FIG. 9. The band structure of Ni along T = (0, 0, 0) and X = (1, 0, 0) and along T and 
L — (0.5, 0.5, ,0.5). The solid curves are the experiment and the dotted curves are the LDA (from 
Martensson and Nilsson 1984). The filled circles are the quasiparticle energies in the GWA. After Aryase- 
tiawan (1992a) and Aryasetiawan and von Barth (1992b). 

As mentioned before, the 6 eV satellite is not reproduced in the GWA. That the 6 eV satellite is 
missing in the GWA can be seen directly in the imaginary part of the self-energy. For the existence 
of a satellite, there should be a strong peak at around 5-6 eV reflecting the presence of a stable 
excitation. But as can be seen in figure ^ such peak is missing. 

3p-resonance photoemission measurements at 67 eV photon energy corresponding to the binding 
energy of the 3p core exhibit an asymmetric (Fano) resonant enhancement of the satellite and the 
main 3d line shows a strong antiresonance (Guillot et al 1977). This is explained as an Auger 
process where a 3p core electron is excited to fill the empty 3d states followed by a super Coster- 
Kronig decay 



3/3d^4s + hLj-^ 3p^3di"4s -^ 3/3^^45 + e 



(148) 



Although there is some indication that the 6 eV structure might arise from single-particle states 
(Kanski, Nilsson, and Larsson 1980) the 3p resonance photoemission provides a strong evidence of 
the many-body character of the 6 eV satellite. The standard explanation for the presence of the 
satellite (Penn 1979, Liebsch 1979, 1981) is that a 3d hole created in a photoemission experiment 
introduces a strong perturbation due to its localized nature, causing another 3d electron to be 
excited to an empty state just above the Fermi level. In atomic picture, the state with two 3d 
holes correspond to the configuration 3cf4:S^ which is separated from the main line configuration 
3(P4:S by more than 6 eV but which is reduced considerably by metallic screening. The two holes 
scatter each other repeatedly forming a "bound state" at 6 eV. In a simple picture, the photon 
energy is used to emit a d electron and to excite another into an empty d state so that the emitted 
electron appears to have a lower binding energy (satellite). The excited electrons mainly come 
from the bottom of the d band since they have the largest mixing with the s-p states and therefore. 
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according to the dipole selection rule, a large transition to the empty d states. This then has been 
argued as the source of band narrowing. The GW calculations, however, show that the largest 
contribution to band narrowing comes mainly from screening rather than two- hole interactions. 



1,0 



0,0 



(a) 


1 




■ 






1 


minority spirt 


r'25 




1 Im 1 1 f 


I, 






^ 


J\ 


, Jl,. 




/Av — " — ^^ 













-10.0 



0.0 

W (a.Li.) 



10.0 




-1(1.1) 



0.0 
(i>(a.u.) 



1 0.0 



FIG. 10. The real and imaginary parts of the correlation part of the self-energy for the minority and 
majority spin state T2S (The lowest 3d state at the F point). The units are in a.u. (Hartree=27.2 eV). 
After Aryasetiawan (1992a). 

The reduction in the exchange splittings is related to the satellite. For simplicity we consider 
a two-level model with fully occupied majority channel and one occupied minority channel. A 
photocmission hole in the majority channel can induce another hole in the minority channel but 
not vice versa since the majority channel has no empty states. Thus, the effects of two-hole 
interactions are larger for the majority than for the minority channel resulting in a reduction in 
the exchange splitting. Calculations based on the Hubbard model within the t-matrix approach 
(Kanamori 1963) confirm this picture (Liebsch 1979, 1981, Penn 1979). These model calculations 
also assign the reduction in the band width as originating from the two-hole interactions which is 
not in complete agreement with the GH^ calculations. These results may be reconciled as follows. 
In the t-matrix calculations, it was found that there was no value of U that gave the correct 
satellite position and the band width. To get the correct band width, the value of U was such 
that the satellite energy became too large (Liebsch 1981). If we assume that the t-matrix theory 
gives the correct physical description for the satellite, the appropriate value of U would not give 
a large reduction in the band width, but this is taken care of by the GWA. Thus one concludes 
that self-energy calculations which include the diagrams of the GWA and t-matrix theory may give 
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both the correct band width and the satelhte structure. Diagrammatic comparison between the 
GWA and the t-matrix theory reveals that hole-hole interactions described by the t-matrix are not 
included in the GWA except to second order only. Direct comparison between G W calculations 
on real systems and Hubbard model calculations is, however, difficult if not impossible. This is 
because the Hubbard U cannot be easily related to the screened interaction in the GW^ calculations. 
An extension of the t-matrix theory including Fadeev's three-body interaction (Fadeev 1963) 
was made by Igarashi (1983, 1985) and by Calandra and Manghi (1992, 1994). The theory has 
been applied to Ni (Igarashi et al 1994, Manghi, Bellini, and Arcangeli 1997) and NiO (Manghi 
and Calandra 1994). 



2. Nickel Oxide 

NiO is a prototype of the Mott-Hubbard insulators. It was pointed out by Mott in the late forties 
that a system with an on-site Coulomb energy larger than the single-particle bandwidth tends to 
become an insulator and that single-particle theory is bound to give a wrong prediction for the 
state of the system (Mott 1949). Indeed, the LDA predicts NiO to be a metal when the calculation 
is performed in a paramagnetic state (Mattheiss 1972a, b). Slater (1974) suggested that a gap could 
be opened by an interplay between antiferromagnetism and crystal-field splittings. A detailed work 
along this direction can be found in a paper by Tcrakura et al (Terakura et al 1984). The LDA does 
produce a gap in an antifcrromagnetic state but of only 0.2 eV in contrast to the experimental gap 
of 4.0 eV (Powell and Spicer 1970, Hiifner et al 1984, Sawatzky and Allen 1984). As expected, the 
free-electron like O p band is well described by the LDA but the magnetic moment is too small (1.0 
/Ub) compared to experiment (1.7-1.9 ^b) (Alperin 1962, Fender et al 1968, Cheetham and Hope 
1983). Clearly there is something seriously wrong with the LDA. A more convincing evidence is 
provided by CoO, where the number of electrons in the paramagnetic structure is odd, making it 
impossible for any single-particle theory to predict CoO as an insulator without doubling the unit 
cell. Experimentally, CoO is an antiferromagnetic insulator so that one might argue like Slater 
that single-particle theory could still give the correct prediction if the calculation is performed 
in an antiferromagnetic structure. One realizes, however, that the difference in magnetic energy 
distinguishing the paramagnetic and antiferromagnetic states is only a fraction of an eV, which 
is much smaller than the band gap. This means that the result of a theoretical calculation for 
the band gap should not depend on whether the calculation is performed in a paramagnetic or an 
antiferromagnetic structure. 

The basic physics of the Mott-Hubbard insulators was explained by Mott several decades ago 
(Mott 1949). From the tight-binding limit, switching on hopping matrix elements causes the 
formation of a band of width w centred around the atomic eigenvalue. The possibility of occupying 
states with lower energy favours electron hopping but it costs a Coulomb energy U for an electron 
to hop from one site to the neighbouring site. If U is larger than w, the gain in kinetic energy is 
overwhelmed by the loss in Coulomb energy and the system prefers to be an insulator with a gap 
approximately given by C/, splitting the lower and upper Hubbard bands. While the Mott picture 
is essentially correct, there are a number of experimental data which cannot be explained. The 
value of [/, for instance, is estimated to be 8-10 eV which is much larger than the experimental gap. 
More recent studies initiated by Fujimori, Minami, and Sugano (Fujimori et al 1984) and Sawatzky 
and Allen (Sawatzky and Allen 1984) based on the cluster approach and Anderson impurity model 
show that the gap in NiO is a charge-transfer gap. If an electron is removed from a Ni site the 
number of holes increases leading to a state with high energy due to an increase in the Coulomb 
interaction among the holes. The hole created on the Ni site may be filled by the transfer of 
an electron from an O site. Although it costs some energy transfer this leads to a state with 
a small binding energy . The states at the top of the valence band therefore have a large O p 
character. The lowest conduction state is of d character as in the Mott picture and the gap is 
therefore formed between the valence O p and conduction Ni d states. Much of the Ni d weight 
goes into a satellite located below the O p, opposite to the Mott and the Slater pictures. This 
model is able to explain experimental data which would otherwise be difficult to explain by the 



47 



Mott and Slater pictures. The most convincing evidence supporting the charge-transfer picture is 
the 2p resonant photocmission experiment in CuO (Tjeng et al 1991) which has a similar electronic 
structure as that of NiO. In this experiment, a 2p core electron is excited and the remaining hole is 
subsequently filled by a valence electron. Since dipole transition matrix element is largest between 
p and d states, resonance in the valence energy region can be identified as the position of d states 
which turns out to be below the O p band, rather than above as in the Mott and Slater pictures. 
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FIG. 11. (a) Comparison between the LDA (solid line) and the experimental band structure of Antifer- 
romagnetic II NiO. From Shen et al (1990, 1991a, 1991b). (b) Comparison between the LDA band structure 
(solid line) and the quasiparticle band structure in the GWA (filled circles) for NiO. After Aryasetiawan 
and Gunnarsson (1995). 



The LDA antiferromagnetic bandstructure is shown in figure |ri|. The highest valence state is 
formed by the majority Cg and minority t2g Ni states and the lowest conduction band is formed 
by the minority eg state. GW calculation for NiO gives a gap of ~5.5 eV (Aryasetiawan and 
Gunnarsson 1995). Starting from the LDA antiferromagnetic bandstructure with a gap of 0.2 
eV, a GW calculation shifts the Ni Cg conduction upwards, increasing the gap to ~1 eV. This 
upward shift of the Cg conduction band leads to a substantial change of the character of the 
wavefunctions, reducing the amount of minority eg character in the occupied states. To include 
a limited self-consistency, we introduce a new non-interacting Hamiltonian H^ with a non-local 



potential Ajcg 



where A is chosen so that the band gap obtained from H" agrees with the 



gap obtained from the previous GW calculation. This i7° is then used to generate a new G*^ 
and a new self-energy in the GWA. This procedure is iterated to self-consistency. The non-local 
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potential modifies the eigenvalues as well as the wave functions used to construct the zeroth order 
Green function G*^. The raising of the unoccupied majority Cg band by the self-energy correction 
reduces the hybridization with the O p band and has the effects of raising the bottom of the O 
p band and pushing down the top of the O p band at the T -point resulting in better agreement 
with photoemission data. In addition, the width of the unoccupied Cg band is reduced. The 
reduction in hybridization also reduces the magnitude of the exchange interaction of the eg band 
with the occupied states which has the consequence of widening the gap. Thus, it is important 
that the wave functions are also modified in the self-consistent procedure. The final position of the 
unoccupied eg band is just below the Ni 4s. As a check, the calculation has also been performed 
in the ferromagnetic state. A gap of ~ 5.2 eV was obtained, close to the antiferromagnetic value. 
In contrast to the Slater model, the gap does not depend on the antiferromagnetic ordering and 
the results correctly predict that NiO remains an insulator above the Neel temperature. The 
miscalculation clearly improves the LDA gap markedly and it is in reasonable agreement with 
the experimental value of 4.0 eV. An estimate of the magnetic moment yields a value of 1.6 /is 
(Aryasetiawan and Gunnarsson 1995) in good agreement with the experimental value of 1.7-1.9 
Mb- 




FIG. 12. The spectral function of antiferromagnetic NiO at k = |(1 0)27r/a within the GWA 
projected on the O p orbitals (solid line) and the Ni d orbitals (The dotted line corresponds to the 
majority spin and the dashed line the minority spin). After Aryasetiawan and Gunnarsson (1995). 

To study the character of the states at the top of the valence band, the spectral function has 
been calculated. Projection of the spectrum into the Ni 3d and O p orbitals shows that there is an 
increase of the O p character at the top of the valence band but the main character is primarily Ni 
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3d. The satellite at -8 eV (Shcn et al 1990, Shen et al 1991a, Shen et al 1991b) is not reproduced 
at the final self-consistent spectrum but a more detailed study of the satellite structure reveals 
that it is rather sensitive to the starting Hamiltonian (Aryasetiawan and Karlsson 1996). Starting 
from the LDA Hamiltonian in fact gives a satellite at about -10 eV but this satellite diminishes 
in intensity as the gap opens up. The origin of this behaviour can be traced back to the presence 
of a plasmon-like peak at low energy which is related to the incorrect LDA bandstructure. As 
the gap opens up, this peak structure becomes broadened and consequently the satellite structure 
diminishes. As in the case of Ni, it appears that the satellite structure is due to short-ranged 
correlations which are not properly taken into account by the RPA. T-matrix approach could be 
appropriate and might remove some Ni d weight from the top of the valence band to the satellite 
region but it is not clear how this could increase the charge transfer from the O p. 



3. Manganese Oxide 

A calculation on MnO based on a simplified GW scheme described in section 4 has been per- 
formed by Massidda et al (1995a). The electronic structure of MnO is the simplest among the 
transition metal oxides and in some respects similar to NiO. As in the case of NiO, the LDA gives 
a too small band gap of 1.0 eV compared with the experimental value of 3.8-4.2 eV. The magnetic 
moment is also somewhat too small (LDA 4.3 /ie, experiment 4.6-4.8 /zb) although the relative 
discrepancy is not as large as in NiO. The larger LDA band gap is due to the fact that in MnO 
the majority spin is fully occupied and the minority spin is empty resulting in a large magnetic 
moment so that the exchange splitting is also large and dominates the ligand-field splitting and 
band broadening due to intersublattice coupling. In NiO the magnetic moment is smaller and the 
exchange splitting is comparable to the ligand-field splitting and band broadening. 

The semiempirical model GFF calculation gives a gap of 4.2 eV which compares well with the 
experimental value of 3.8-4.2 eV. The LDA magnetic moment is also improved to 4.52 jx-q- There is 
an increase of O p character and a decrease of Ni 3d character at the top of the valence band, which 
are percentagewise large but small in absolute term, so that the main character is still primarily 
Ni 3d. The results are qualitatively similar to the full Giy calculations on NiO described above. 
Calculation on CaCu02 using the model GW scheme was also done recently (Massidda et al 1997). 



4- 3d semicore states 

It is well-known that the LDA eigenvalues for localized states are usually too high compared to 
experiment. The discrepancies can be several eV. For example, the Zn semicore 3d states in ZnSe 
are too high by 2.5 eV, the Ga semicore 3d states in GaAs by 4 eV and the Ge 3d semicore states 
by as much as 5 eV. In free atoms these deviations are even larger and in, for instance, a free Zn 
atom the 3d eigenvalue is about 6.5 eV too high. This raises interesting questions about whether 
or not GWA can describe these shallow core levels and why the LDA eigenvalues are substantially 
worse for free atoms than for solids, although the semicore states are almost completely localized. 

For deep core levels, an important contribution to the GWA self-energy comes from the polariza- 
tion (relaxation) of the more weakly bound electrons (Hedin and Lundqvist 1969, Lundqvist 1969). 
This relaxation is a classical electrostatic effect, which can be described in ASCF calculations, per- 
forming ground-state calculations with and without the core hole. Explicit ASCF calculations for 
free atoms were performed by Hedin and Johansson (1969), who obtained quite accurate results. 
This suggests that GWA may be rather accurate for such deep core levels, since it includes similar 
physics as the ASCF calculations. 
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FIG. 13. Partial densities of spin-up states of antiferromagnetic MnO. Top panel: the total O 2p 
projection into the O spheres. Second and third panels: d projection into the two inequivalent Mn spheres. 
Solid (dashed) lines correspond to the model GW (spin polarized LDA) calculations. Bottom panel, up- 
per curves: Inverse photoemission spectrum (from van EIp et al 1991) and the difference between the 
on- and off-resonance photoemission spectra (from Lad and Henrich 1988), representing the Mn contribu- 
tion. Lower curves: GW d-projected density of states into both Mn spheres. Dotted lines illustrate the 
integrated-intensity background. After Massidda et al. 

It is not a priori clear that GWA can describe semicore states. Unlike the deep core states, the 
shallow core states are screened by states which are almost degenerate with the hole, which may 
introduce new effects. Calculations for semicore states in ZnSe, GaAs and Ge (Aryasetiawan and 
Gunnarsson 1996) showed, however, that GWA improves the LDA eigenvalues significantly also 
for these core states, leaving a discrepancy of only 0.5-1.0 eV to experiment (too small binding 
energy). 

The errors in the LDA eigenvalues are to a large extent due the unphysical interaction of an 
electron with itself. This interaction is only incompletely cancelled by the exchange correlation 
potential (Gunnarsson, Lundqvist, and Wilkins 1974). This leads to an ionization energy which is 
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too small. In addition relaxation effects are neglected, which tends to overestimate the ionization 
energy and therefore partly cancels the error from the self-interaction. As discussed above, essential 
relaxation effects are included in ASCF calculations. In addition, the LDA total energies are rather 
accurate and the self-interaction errors are rather small. This is due to an exact sum-rule satisfied 
by the LDA (Gunnarsson and Lundqvist, 1976). While this sum-rule is very important for total 
energy calculations, it does not necessarily imply a good exchange-correlation potential or accurate 
eigenvalues. It is therefore crucial that ASCF calculations involve total energy differences rather 
than eigenvalues. 

To understand the different accuracy of the LDA eigenvalues for a solid and a free atom, Aryase- 
tiawan and Gunnarsson (1996) performed transition state calculations in the Slater transition state 
approach (Slater 1974). These calculations showed that in the solid the creation of a core hole leads 
to a substantial charge transfer to the site where the hole was created (Lang and Williams 1977, 
Zunger and Lindefelt 1983). Due to this charge transfer the relaxation energy is larger than for a 
free atom. The cancellation of the error in the eigenvalues due to the self-interaction and the error 
due to the neglect of relaxation effects therefore becomes more complete in the solid than in a free 
atom. This explains why the errors in the eigenvalues are much smaller in the solid (Aryasetiawan 
and Gunnarsson 1996). 

The ASCF gave ionization energies for the 3d semicore states which were about 1 eV too large. 
The discrepancy between ASCF and GWA may be due to an error in the LDA, which tends to 
overestimate the exchange-correlation energy between a 3d electron and the 3s-3p core (Gunnarsson 
and Jones 1980) which is about 1 eV in the series K-Cu (Harris and Jones 1978). It could also 
be due to the difference in the RPA screening and the static LDA screening. Transition state 
calculations for a number of Zn compounds have also been performed by Zhang, Wei, and Zunger 
(1995). We observe, however, that the Slater transition state approach requires that the semicore 
states are sufficiently localized to form a bound state in the transition state calculations. The 
bound state has no dispersion, in contrast to the GW^ results which show the full band structure. 
Thus, if a bound state is not formed in the transition state calculation, the result would be identical 
to that of the LDA. The problem can be illustrated for Cu metal where the 3d band is 0.5 eV too 
high. The transition state would give a single number rather than a band assuming that a bound 
state is formed in the first place. A G M^ calculation on the other hand lowers the position of the 
band by 0.2-0.3 eV while maintaining the LDA band structure (Aryasetiawan and Gunnarsson 
1997). 

A model GW^ calculation for ZnO by Massidda et al (1995b) also improved the LDA result from 
-5.4 eV to -6.4 eV but a significant error remains when compared to experimental results -8.6 eV 
and -7.5 eV. 



D. Surfaces 

1. Jellium surface 

GW calculations for a model jellium- vacuum interface have been done in some details. The 
problem of interest here is how the surface barrier looks like. According to classical electrostatic, 
the image potential seen by an electron in the vacuum far from the surface behaves like 



1 

4 (z - zo) 



^image — ^ ^ ^^ (^i4yj 



where z is the coordinate normal to the surface and zq is the position of the effective image plane. 
In the LDA, the image potential is known to decay exponentially (Lang and Kohn 1973). This 
unphysical behaviour is due to the fact that in the LDA the exchange-correlation hole is determined 
by the local density and it does not feel the surface directly. The exchange-correlation hole obeys 
a sum-rule that it must integrate to one but since the density is small outside the surface it means 
that the hole becomes very extended. In fact about half the hole resides far inside the surface. 
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Due to the very extended structure of the hole, the resulting V^"^ decays exponentially outside 
the surface. The situation is similar to the atomic case where the LDA V^"^ decays exponentially 
instead of decaying as l/r as in the exact V'^'^. (Almbladh and von Barth 1985a). Similarly, the 
exact DFT exchange-correlation potential should be capable of reproducing the correct behaviour 
of the image potential which is important in many applications. For example, the LDA potential 
cannot produce the Rydberg series of image states. Binding energies and life-time of surface 
states bound by image potential depend crucially on the 1/z dependence as do tunneling currents 
in the scanning-tunneling microscope and positions of alkali ions adsorbed on metal substrates. 
Interpretation of inverse photoemission data relies on the existence of the image tail of the surface 
barrier. 

By calculating the self-energy for the surface, a V'^'^ can be derived simply by taking the trace 
of the following Dyson equation (Sham and Schliiter 1983,1985, van Leeuwen 1996) 

G ^ G°^ + G°^ (E - V"") G (150) 

and noting that the exact density is equal to the exact DP density (i.e. the trace of G is equal to 
the trace of G^^). This yields 

J d^nV^^n) I dwG'"' {r,r„u;)G{n,r,u;) 

= j d\^ j d\2 j duG""^ {V,V^,UJ)^{VI,V2,L0)G{V2,V,U) (151) 

The exchange-correlation potential obtained from the above equation using the GFF self-energy 
exhibits the correct asymptotic behaviour as can be seen in figure WA It can also be seen from the 
figure that it is the correlation potential originating from the Coulomb correlation that determines 
the asymptotic behaviour of the image potential. The exchange potential, on the other hand, is 
numerically found to decay as v^ ~ —a/z^ and to contribute significantly to the determination of 
zq. The position of the effective image plane zq deduced from V^'^ is therefore different from the 
one for a classical test charge. Thus, for rg — 2.07 the value of zq deduced from the image tail of 
y™ is zq = 0.72 ± 0.1 a.u. measured from the jellium edge while the value obtained from a linear 
response calculation is zq = 1.49 a.u. (Eguiluz et al 1992). 

The Kohn-Sham eigenvalues calculated from the V^^ deduced from Sgvk turn out to very close 
to the quasiparticle energies obtained from the Dyson equation (Deisz, Eguiluz, and Hanke 1993). 
The difference Eqp — Eks — 0-02 eV for q\\ = which is much smaller than the binding energy 
of the state - about 0.5 eV. Furthermore, the Kohn-Sham and quasiparticle eigenfunctions are 
practically identical (the overlap is > 0.999 ). This does not, however, imply that the physics of 
the self-energy at the surface can be completely described by a local potential. The imaginary 
part of the self-energy associated with damping is intrinsically non-local and energy dependent. It 
cannot be mimicked by a local complex potential. As the electron moves away from the surface, 
the maximum of Im S remains at the surface, reflecting a high degree of non-locality. The range of 
this non- locality is comparable with the decay length of the electron density in the vacuum (Deisz, 
Eguiluz, and Hanke 1993). It was also found that Im S deviates from the quadratic behaviour 
~ {E — Ep)^ when the electron is in the vacuum outside the surface and the departures grow as the 
electron-surface separation increases. The effect is attributed to the suppression of one-electron 
decay channels near the Fermi level (Deisz and Eguiluz 1997). 
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FIG. 14. V^[z) at a jellium surface for r^ 
using the GWA for E, ai 
dashed curve is the image potential l/""(z) 



3.93 (Af = 12.9 a.u.). The sohd curve is the solution 
to equation (151) using the GWA for E, and the dotted curve is the corresponding LDA potential. The 

-l/4(z - zo). After Equiluz et al (1992). 



2. Si(l 0) surface 

There have also been studies for more detailed microscopic models of surfaces. In particular, 
the Si(lOO) surface has been extensively studied experimentally and theoretically, because of the 
technological importance of Si. At room temperature there is a 2 x 1 reconstruction which at low 
temperatures goes over in a c(4 x 2) reconstruction (see, e.g., Johansson et al. for references to 
experiment). The surface atoms are believed to form buckled dimers, where one atom moves out 
of the surface and the other into the surface. The room temperature photoemission spectrum has 
been measured by Johansson et al. (1990). 

The electronic structure has been calculated in the GWA by Northrup (1993) (c(4 x 2)), by 
Kress, Fiedler, and Bechstedt 1994 (2 x 1) and by Rohlfing et al. (1995b)(2 x 1). In figure^ we 
show the results of Rohlfing et al. (1995b). These calculations were performed using a Gaussian 
basis set and a super cell containing eight Si layers. The figure shows the LDA bulk bands (dashed) 
and GW surface states. The GW surface band agrees rather well with an experimental band, while 
another experimental band has no correspondence in the theoretical calculation. It is interesting 
that the calculation of Northrup, using the low temperature structure, produced two bands in 
close agreement with experiment. Actually, the experimental samples may have some domains 
with this c(4 x 2) structure (Johansson et. al 1990). The GW calculation shifts the conduction 
band by about +0.50-+0.65 eV relative to the valence band (Rohlfing et al. 1995b). According to 
optical experiments, the bottom of the valence band at F is 1.1 eV above the top of the valence 
band in good agreement with the theoretical result 0.95 eV. For the indirect band gap there are 
experimental estimates in the range 0.44 eV to 0.9 eV compared with the theoretical result 0.7 eV 
(Rohlfing et al. 1995b). 
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FIG. 15. Calculated dangling-bond bands. Solid lines, GWA energies, dashed lines LDA energies. The 
experimental results are shown by diamonds (Uhrberg et al 1981) and circles (filled and open) (Johansson 
et al 1990). (After Rohlfing et al. 1995b). 

In this context we also mention that there has been work within the GWA for interfaces, e.g., 
for calculating Schottky barriers. This work will not be discussed further here, but we refer the 
reader to Charlesworth et al. (1993) and Godby and Sham (1994) for an example of such work 
and further references. For other work on semiconductor surfaces we refer to Hybertsen and Louie 
(1988b) and Bechstedt and Del Sole (1990). Rohlfing, Kriiger, and PoUmann (1996) performed 
calculations on clean, H and S terminated Ge(OOl) surface and Zhang et al (1988) have calculated 
the valence band off-set of AlAs-GaAs(OOl). 



E. Clusters 

Electronic excitations and optical spectra in clusters have been studied mainly within the con- 
figuration interaction approach (Bonacic-Koutecky et al 1992, Bonacic-Koutecky, Fantucci, and 
Koutecky 1990a, b). While it gives accurate results, its applications are limited to systems contain- 
ing a small number of electrons, typically less than 10. GWA provides an alternative for calculating 
excitation properties in clusters with relatively large number of electrons which cannot be handled 
by the CI method. 

GW calculations have been performed for a jellium-sphere model for alkali metals (Saito et al 
1989) and more recently for the real cluster Na4 (Onida et al 1995). In alkali-metal clusters, it is 
known that the ionization energies calculated within the LDA are too low compared to experiment 
(Ishii, Ohnishi, and Sugano 1986, Saito and Cohen 1988) and the discrepancy becomes worse 
the smaller the cluster. The size dependence of the ionization energy in the LDA is too weak. 
This discrepancy is attributed to self-interaction which is not taken into account properly in LDA 
eigenvalues. That the discrepancy gets worse for smaller clusters is intuitively clear since the larger 
the clusters the more they resemble electron gas on which the LDA is based. On the other hand, 
the LDA gives the correct sequence of eigenvalues for the valences states in the jellium-sphere 
model which is Is, Ip, Id, 2s, 1/, 2p, etc. (Saito et al 1989). 

In the jellium-sphere model, the ionic charges are smeared to form a sphere of a uniform positive 
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background. GW calculation for this model corresponding to Na2o lowers the occupied LDA 
eigenvalues and increases the unoccupied ones, thus increasing the energy gap similar to the results 
for bulk semiconductors and insulators. The ionization energy as a function of cluster size is also 
in better agreement with experiment although the absolute values are somewhat too large (Saito et 
al 1989). This could be due to the absence of core polarization in the jellium-sphere model which 
gives a positive contribution to the self-energy. The GW results are actually rather close to the 
HF results, the reason being that screening due to long-range correlations is much less important 
in a small finite system than in a solid (Saito et al 1989). 

G W calculation for real Na4 yields similar results (Onida et al 1995). The unoccupied states 
are raised by between 0.75 and 0.90 eV while the highest occupied molecular orbital (HOMO) and 
the lowest occupied molecular orbital (LUMO) state are lowered by 1.55 and 1.40 eV respectively 
giving a HOMO-LUMO gap of 3.0 eV compared with the LDA value of 0.55 eV. Again the GW 
result is close to the HF value of 3.4 eV for the reason discussed above. Indeed calculation of the 
static dielectric function for this cluster shows a metallic behaviour for small distances but it drops 
to less than unity starting at about 6 a.u.. For comparison, the value of the dielectric function 
for metallic sodium at 5 a.u. is about 50 whereas in the cluster it is 1.2. This phenomenon of 
antiscreening is typical of a small system. If a positive test charge is introduced at the centre of the 
cluster, say, electrons will surround the test charge and the surface of the cluster will therefore be 
positively charged. Thus as the distance from the centre increases, the screening quickly vanishes 
and becomes negative (Onida et al 1995). A similar behaviour is also observed in the Cgo molecules 
(Gunnarsson, Rainer, and Zwicknagl 1992). 
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FIG. 16. Calculated absorption spectra of Na4 including self-energy and excitonic effects (dotted 
line) in arbitrary units and using a Gaussian broadening of 0.06 eV. The solid line is the experimental 
photodepletion cross sections from Wang et al (1990). The vertical bars show the unbroadened spectrum 
and the inset shows the LDA results. After Onida et al (1995). 



F. Fullerenes 

Solid Ceo (fuUerite) is a molecular solid, where the orbitals of a free Cgo molecule essentially 
keep their character and the hopping between the molecules only leads to a small broadening of 
the discrete molecular states. The band structure of Cgo is shown in figure O. In undoped Cgo 
the Hu band is occupied and the Tiu band is empty. The LDA (left part of the figure) band 
gap is about 1 eV (Erwin and Pickett 1991, TrouUier and Martins 1992, Satpathy et al. 1992), 
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which is substantially smaller than the experimental value, 2.3 eV (Lof et al. 1992), obtained from 
photocmission and inverse photoemission. 
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FIG. 17. LDA (a) and GW (b) band structures for solid Ceo in the Fm3 structure 
(After Shirley and Louie 1993). 

Shirley and Louie (1993) have performed a GW calculation for solid Ceo on a fee lattice with all 
molecules having the same orientation {Fm3 structure). This differs somewhat from the experi- 
mental T = structure Pa3, where the molecules take four different orientations. This difference 
should not be important for the present discussion. Shirley and Louie (1993) used the Levine and 
Louie (1982) model for the static dielectric function together with a plasmon pole approximation. 
Their results are shown in figure |l^. The band gap is increased to 2.15 eV in good agreement 
with the experimental results. The band widths were also increased by about 30 %. There are 
no reliable experimental results for the dispersional band widths. It was concluded that undoped 
solid Ceo is a standard band insulator with dispersive bands. 

The calculation by Shirley and Louie was performed for the undoped Cgo solid. In A3 Ceo (A= 
K, Rb) the Ti„ band is half-filled. The Hubbard U ~ 1.5 eV (Lof et al. 1992), describing the 
interaction between two electrons on the same molecule, is large compared with the width (W ~ 0.6 
eV) of the partly filled Tiu band. One may therefore expect strong correlation effects for these 
systems. Although it is not clear if the GWA can describe such strong correlation effects, it is 
interesting to perform such a calculation. 

To screen the Hubbard U, Gunnarsson (1997) considered a model dielectric function 



e(q,w) =eo- 



w, 



Op 



(152) 



which describes the coupling to a plasmon at uj — uJop/\/m*eo. This plasmon is due to t he 
oscillations of the three electrons donated by the alkali atoms to the Ti^ band. The model (152) is 
therefore complimentary to the calculation of Shirley and Louie (1993), since this Tiu plasmon does 
not exi st fo r the undoped system, while the physics considered by Shirley and Louie is neglected 
in Eq. ( |152| ). This model of the dielectric function was combined with a tight-binding (TB) model 
(Satpathy et al. 1992) for the band structure, which reproduces the LDA Tiu band well. 

First a Hartree-Fock (HF) calculation was performed for this model. It was found that in HF 
the width of the Ti„ band is increased by about 75 %. Next a GW calculation was performed. 
Including the coupling to the Tiu plasmons was found to reduce the Ti„ width to a value 35 % 
smaller than the original TB width. Actually for a model of this type one can show under rather 
general assumptions that the band width is always reduced if there are no other bands above or 
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below the band considered (Gunnarsson 1997). The density of the Ti„ electrons is very small 
and corresponds to the electron gas parameter r^ ~ 7. Thus the density is substantially smaller 
than for the free-electron like metals. It is then not too surprising that the quasi-particle weight 
Z ~ 0.4 — 0.5 is also smaller than for these metals. This means that much of the spectral weight 
appears in satellites. If the narrowing of the Ti„ band due to the coupling to the Tiu plasmon is 
combined with the broadening found by Shirley and Louie due to other couplings, the net result 
is a small change of the band width. 



VI. SELF-CONSISTENCY 

The set of Hcdin's equations (^Q,^,^6|) , in the original formulation of the self-energy expansion 
in powers of the screened interaction W constitutes a self-consistent cycle (Hedin 1965a, Hedin and 
Lundqvist 1969). Within the GWA, starting from a (usually) non-interacting Green function Gq 
one calculates the polarization function Pq = —iGoGa and the corresponding screened interaction 
Wq. The self-energy is then obtained from Sq — iGoWo. In most GW calculations that have 
been performed so far, Eg is taken to be the final self-energy. The interacting Green function 
G obtained from the Dyson equation G = Go + GoTiqG is, however, not necessarily the same 
as Go- To achieve self-consistency, the Green function obtained from the Dyson equation should 
be used to form a new polarization function P — ~iGG, a new screened interaction W, and a 
new self-energy S which in turns yields a new Green function through the Dyson equation. This 
process is continued until G obtained from the Dyson equation is the same as G used to calculate 
the self-energy. Self-consistency is evidently an important issue since it guarantees that the final 
results are independent of the starting Green function. Moreover, according to the Baym-Kadanoff 
theory (1961), a self-consistent GW scheme ensures conservation of particle number and energy 
when the system is subjected to an external perturbation. Conservation of particle number means 
that the continuity equation 

-dtn{r,t)=V-Jir,t) (153) 

is satisfied. Conservation of energy means that the energy change when an external potential 
is applied to the system is equal to the work done by the system against the external potential 
when calculated using the self-consistent G (Baym and Kadanoff 1961, Baym 1962). Moreover, 
self-consistency ensures that 

1 f 
N = -tr dujlmG{uj) (154) 

T^ J-oo 

gives the correct total number of particles, when n and j is obtained from the self-consistent Green 
function. The first self-consistent GW calculation was probably by de Groot, Bobbert, and van 
Haeringcn (1995) for a model quasi-one dimensional semiconducting wire. The relevance of this 
model to real solids is, however, unclear. 

Five aspects may be distinguished in relation to self-consistency (von Barth and Holm 1996): 

1. Modification of quasiparticle wavefunctions. 

2. Shift of quasiparticle energies. 

3. Modification of quasiparticle weights (Z factors). 

4. Modification of quasiparticle life-times. 

5. Modification of the screening properties of the system. 

These aspects were recently studied in detail by von Barth and Holm (1996) and by Shirley 
(1996) for the electron gas. Naturally, the first aspect cannot be addressed for the electron gas 
since the quasiparticle wavefunctions remain plane waves. The results of these studies are 
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• The band width is increased from its non-self-consistent value, worsening the agreement with 
experiment. 

• The weight of the quasiparticles is increased, reducing the weight in the plasmon satellite. 

• The quasiparticles are narrowed, increasing their life-time. 

• The plasmon satellite is broadened and shifted towards the Fermi level. In fact, it almost 
disappears at full self-consistency. 
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FIG. 18. The quasiparticle dispersion for the electron gas for r^ = 2 and 4 from partial self-consistent 
GWo calculations (Wo is held fixed at the RPA whereas G is allowed to vary to self-consistency) compared 
with the free-electron dispersion. The largest change in the band width occurs for r^ = 4. After von Barth 
and Holm (1996). 

As discussed below, the main effects of self-consistency are caused by allowing the quasiparticle 
weight to vary (von Barth and Holm 1996). These results are true for the case when the screened 
interaction W is fixed at the RPA level (calculated using the non- interacting Gq) and only the 
Green function is allowed to vary to self-consistency and also for the case when both G and W are 
allowed to vary (full self-consistency case). The increase in the band width is disturbing and can 
be understood as follows. Let us consider the first case with fixed W = Wq for simplicity. First 
we note that the GW result for the band width after one iteration is close to the free electron one. 
This means that there is almost a complete cancellation between exchange and correlation. After 
one iteration the quasiparticle weight is reduced to typically 0.7 and the rest of the weight goes 
to the plasmon satellite. Correspondingly, when the new G is used to calculate a new Im S"^, its 
weight is transferred away from low energy to high energy, due to the sum rule (von Barth and 
Holm 1996) 



dw|lmE^(k,tj) 



E 

q 



dw|ImWo (q, w) 



(155) 



which shows that the left-hand-side is a constant depending only on the prechosen Wo but inde- 
pendent of k and self-consistency. For a state at the Fermi level, this has little effect since Im S'^ 
has almost equal weights for the hole {oj < /i) and the particle part {uj > /i) which cancel each 
other when calculating ReS'^, as can be seen in equation (pq). But for the state at the bottom of 
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the valence band, Im S"^ has most of its weight in the hole part so that the shifting of the weight in 
Im Yi'^ to higher energy causes Re E"^ to be less positive than its non-self-consistent value. A similar 
effect is found for the exchange part which becomes less negative but because the bare Coulomb 
interaction has no frequency dependence, the renormalization factor has a smaller effect on S^ so 
that the reduction in S^ is less than the reduction in S"^. The net effect is then an increase in the 
band width. The shifting of the weight in IniS'^ to higher energy has immediate consequences of 
increasing the life-time and the renormalization weight Z (through a decrease in | 9ReE'^/9a>|) of 
the quasiparticles and of broadening the plasmon satellite, compared to the results of one iteration 
(von Barth and Holm 1996). 

When W is allowed to vary (full self-consistency) the results become even worse: the band width 
becomes even wider and the plasmon satellite becomes broad and featureless, in contradiction to 
experiment. The quasiparticle weight is increased further. These results can be explained by the 
disappearance of a well-defined plasmon excitation in W. The quantity P — —iGG no longer has 
a physical meaning of a response function, rather it is an auxiliary quantity needed to construct 
W. Indeed, it does not satisfy the usual /— sum rule. The equations Ree(q, uop) — Ime(q, uup) — 
determining the plasmon energy are not satisfied any more since the Green function now always 
has weight around txi = ujp. This has the effect of transferring even more weight in ImS^ to higher 
energy with the consequences discussed in the previous paragraph. Shirley (1996) included the 
second order self-energy diagram (vertex correction) and found that it cancelled the effects of 
self-consistency to some degree. 
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FIG, 19, The spectral function from a partial self-consistent GWo calculation (see previous figure) 
compared to that of the first iteration for k = Q.^kp and r^ = 4. The self-consistent quasiparticle energy 
is lowered (band broadening) compared to the non-self-consistent one whereas the satellite position is 
somewhat improved. After von Barth and Holm (1996). 

A very interesting outcome of the full self-consistent calculation (Holm and von Barth 1997) is 
that the total energy calculated with the Galitskii-Migdal formula (1958) turns out to be strikingly 
close to the total energy calculated from a much more elaborate quantum Monte Carlo technique 
(Ceperley and Alder 1980). For r^ = 2 and 4 quantum Monte-Carlo gives 0.004 and -0.155 Rydberg 
respectively while the self-consistent GW gives 0.005 and -0.156 Rydberg. This unexpected result 
may be related to the fact that the self-consistent GW scheme is energy conserving and it is partly 
explained by consideration of the Luttinger-Ward energy functional (1960) which is variational 
with respect to G and its minimum is equal to the value obtained from the Galitskii-Migdal 
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formula. What is not clear is why the first order energy diagram (giving the GW self-energy 
upon taking functional derivative with respect to G) appears to represent a very good energy 
functional. Furthermore, the chemical potential calculated from n = dE/dN is in agreement with 
the value obtained from ^ = kp/2 + E{kF,^J') and the particle density n = 2^,^./^^ du;A(k,u;) 
yields n = fc|./(37r^) , i.e. particle number is conserved, as proven by Baym (1962). It has also 
been proven that with a fixed W — Wq particle number is also conserved (Holm and von Barth 
1997). 

The conclusion is that it is not a good idea to perform fully self-consistent GW calculations for 
quasiparticle energies. If full self-consistency is not introduced, an important question is how to 
choose Hq determining Gq. Farid (1997) argued that Hq should be chosen so that its ground-state 
density is the same as the one resulting from the GW calculation. It could be more favourable 
to perform partially self-consistent calculations by, for instance, fixing W at the RPA level and 
modifying the quasiparticle energies in Go, or by choosing a single-particle Hamiltonian such that 
the resulting quasiparticle energies obtained from a GW calculation are the same as the those 
of the single-particle Hamiltonian. In any case, efforts should be directed towards finding vertex 
corrections (beyond GW ). 



VII. VERTEX CORRECTIONS 

By vertex corrections we mean corrections to the self-energy beyond the GW-KPA approximation 
and corrections to the response function beyond the RPA. Vertex corrections to the RPA response 
function naturally involve interaction between screening electrons not taken into account in the 
RPA. This interaction must include the effect of exchange and correlation, which makes the inter- 
action depend upon the relative spin states of the electrons. Vertex corrections to the self-energy 
on the other hand involve interaction of the hole with its surroundings not taken into account in 
the GWA. A consistent treatment should include vertex corrections in both the response function 
and the self-energy (Ward 1950). This preserves conservation laws and is a natural outcome of 
the perturbation expansion of the self-energy in the functional derivative technique of Baym and 
Kadanoff (1961) and Baym (1962). For example, calculations of the optical spectra of Si indicate 
a cancellation between vertex corrections in the response function and the renormalization of the 
quasiparticle energies (Bechstedt et al 1997, Del Sole and Girlanda 1996). Mahan and Sernelius 
(1989) found a similar effect in the calculations of band widths of the electron gas. Conserving 
approximations including vertex corrections with exchange effects only have been considered by 
Hong and Mahan (1994) for the electron gas. 

In the GWA, the screened interaction is calculated within the RPA which takes into account 
primarily long-range correlation which gives rise to collective excitations (plasmons). The pho- 
toemission hole is coupled to one plasmon only in the GWA. Vertex corrections may be loosely 
divided into the short-range and long-range parts. Short-range vertex corrections improve the 
description of the quasiparticles and low-energy satellites whereas long-range vertex corrections 
improve description of high-energy satellites (plasmons). 

The RPA pair distribution function, which is the probability of finding another electron at a 
certain distance from a given electron, is negative for small distances, which is unphysical (Pines 
1961). Exchange and correlations should therefore increase the probability of finding another 
electron at small distances. The same applies to holes. For a test charge, the screening becomes 
more effective when the effects of exchange and correlations are taken into account. But for an 
electron, the RPA screening is reduced when the effects of exchange and correlations between the 
screened electron and the surrounding electrons are taken into account. The physical interpretation 
of vertex corrections is as follows: an electron pushes away other electrons in its vicinity creating a 
screening hole around the electron. Taking into account exchange and correlations of the screening 
holes increases the screening since the probability of the holes getting closer together is increased 
leading to more screening holes. If we now take into account the effects of exchange and correlations 
between the screened electron and the other electrons, then screening will be reduced because the 
probability of finding electrons at small distances to the screened electron is increased, leading to 
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stronger effective interactions between the electrons. The net effect is that the RPA screening is 
reduced. Thus, vertex corrections for electrons will in general reduce screening. Vertex corrections 
in the response function can also take the form of interaction between electrons and holes in 
electron- hole pairs created by the perturbation in the system due to the presence of a photoemission 
hole. This interaction may actually create an additional bound state (exciton) much lower in energy, 
of the order of a few eV, than the plasmon energy. 

Vertex corrections to the self-energy are particularly important for systems with localized states, 
such as those containing 3d and 4f orbitals. This is because the electron-hole pairs created in 
the screening process can be rather localized and therefore interact strongly with the localized 
photoemission hole. This interaction can significantly modify the quasihole energy and its weight 
as well as creating new collective excitations appearing as low energy satellites. This type of 
vertex corrections is short range in nature. Long-range vertex corrections modify the structure of 
the plasmon satellites and may create multiple plasmons as observed in the alkalis. In the GWA, 
the photoemission hole only couples to one plasmon, resulting in one-plasmon satellite. The hole, 
however, can in general interact with several plasmons during its propagation, producing multiple 
plasmon satellites. 



A. Direct evaluation of the second-order self-energy 

Formally, the GWA is the first-order term in the expansion of the self-energy in the screened 
interaction W. It is then thought that the simplest vertex correction is the second-order self-energy. 
This procedure, however, warrants some precautions. First, the physical meaning of the second- 
order term is not clear. Second, this second-order term when evaluated with a frequency-dependent 
interaction can give a self-energy with wrong analytic properties which result in a negative density 
of states, as shown in the electron-gas case (Minnhagen 1974). 

This type of vertex correction was calculated for the band gap of Si both with a bare and a 
screened interaction (Baling and van Haeringen 1989, Daling et al 1991, Bobbert and van Haeringen 
1994). The second-order self-energy with the bare Coulomb interaction was found to be factors 
18 and 37 smaller than the first-order self-energy for Fis and F25 states, which correspond to a 
correction of ~1 % with respect to the Hartree-Fock direct gap. A calculation using a screened 
interaction yields a correction to the gap of the order of ~4 %. Since this calculation was performed 
for the gap states only, it is not clear if the density of states became negative for some energies. 
Both calculation indicate in any case that the second-order term is small. This is encouraging since 
it suggests that higher-order terms are probably small too since the GWA already gives results in 
agreement with experiment. An interesting result is that the second-order vertex correction does 
not shift the absolute position of the LDA valence-band maximum which in the approximation 
used is too low by 0.5 eV. It is speculated that higher order vertex corrections could account for 
the required shift (Bobbert and van Haeringen 1994). 



B. Vertex corrections based on the LDA exchange-correlation potential 

The set of equations (^,^,^,^2) derived by Hedin provides a systematic way of developing 
a perturbation series for the self-energy in powers of the screened interaction W. In the original 
derivation, the zeroth-order Green function is taken to be the Hartree one. However, the band 
structure and wavefunctions in the Hartree theory are less accurate compared with those of the 
LDA, making Hartree Green function a poor starting point. In practical calculations, one uses a 
Green function constructed from the LDA band structure. The response function is also calculated 
using the LDA Green function. 

Consistent vertex corrections were derived by Rice (1965) in relation to the dielectric function 
as discussed in the alkali section. A similar approach was also made Ting, Lee, and Quinn (1975) 
and Mahan (1994). Alternatively, consistent vertex corrections can also be derived from the the 
set of equations ( p3| , p4| ^,46) where one regards V^^^ as a self-energy correction to the Hartree 
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approximation, albeit an ad hoc one (not based on a systematic diagrammatic expansion). Based 
on this starting point, the vertex function A can be easily evaluated yielding (Hybertsen and Louie 
1986, Del Sole, Reining, and Godby 1994) 



A(l, 2, 3) = (5(1, 2)5 (1,3) ~i d (5, 6, 7) K'''' (1, 5) G (5, 6) G (7, 5) A (6, 7, 3) (156) 

where 

SV" (1) 
K-(l,5) = ^ (157) 

remembering that p (1) = —iG (1, 1"*"). The new self-energy with vertex corrections has the form 
of the GWA but with a new screened potential 

W^v[l-P^'{v + if^'^)] "^ (158) 

corresponding to a dielectric function 

?=l-pO(v + if'='=) (159) 

The same result was already discussed in section V.A where the self-energy was expressed as a 



corrected dielectric function e (equation (136)) and a vertex correction F (equation (137)). The 
product e~"'^r is equal to e^^ with e given above. This dielectric function may also be derived 
straightforwardly from time-dependent LDA and may therefore be interpreted as the dielectric 
function that screens the external potential felt by an electron, as opposed to a test charge (Hy- 
bertsen and Louie 1986). This distinguishes itself from the RPA (time-dependent Hartree) dielec- 
tric function in that the induced charge does not only generate the Hartree potential but also an 
exchange-correlation potential. It is clear that one could start with a different local zeroth-order 
self-energy other than V^da ^^'^ arrive at a similar formula. Note that in the formula for the 
self-energy in equation (Esh, the vertex function A enters both W (through P) and S. A problem 
with starting with a local zeroth-order self-energy is that the self-energy with the vertex corrections 
becomes asymmetric in r and r'. 

Application of such scheme to the electron gas gives small changes compared with the original 
GWA (Mahan and Sernelius 1989). Similarly, for Si it yields practically the same gap (0.70 eV) 
and valence band width (11.4 eV) as those of the standard GW calculation (Del Sole, Reining, 
and Godby 1994). However, the absolute position of the top of the valence band is shifted 10 meV 
upwards by the GW-I- vertex corrections and 400 meV downwards by the standard GFK calculation. 
This could improve the calculations of band off-sets at interfaces. A calculation with the vertex 
corrections in the response function R only but not in E, that is with 

W = V \l + V [l - P° {v + K""")]'^ P°] (160) 



R^ [1 - po (t- + X'^^)] 'po, (161) 

gives a smaller gap (0.57 eV) and a smaller band width (10.9 eV) (Del Sole, Reining, and Godby 
1994). The latter result is similar to the case of the alkalis (Northrup, Hybertsen, and Louie 1987, 
Surh, Northrup, and Louie 1988). The top of the valence band is shifted downward by 0.4 eV 
like in the standard GWA. The result for the band width should be taken with caution because 
the plasmon-pole approximation used in the calculations. The worse results obtained by including 
vertex corrections in R only, is consistent with the violation of conservation laws. 
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C. The cumulant expansion 

A diagrammatic approach for including vertex corrections is provided by the cumulant expansion 
method. One of the first applications of the cumulant expansion method was in studying the X-ray 
spectra of core-electrons in metals (Nozieres and de Dominicis 1969, Langreth 1970). Later on the 
method was extended to valence states by Bergersen, Kus, and Blomberg (1973) and Hedin (1980). 
The core-electron problem is modelled by a Hamiltonian consisting of a core electron interacting 
with a set of plasmons: 

H = ec^c + Y, ^q&^^ + Y. ^^^5q (^ + ^q) (162) 

q q 

where c is the annihilation operator for the core electron with energy e, 5j_ is the creation operator 
for a plasmon of wave vector q and energy Wq. and the last term is the coupling of the core electron 
to the plasmon field. The Hamiltonian can be solved exactly and it can be shown that the cumulant 
expansion also gives the exact solution (Langreth 1970): 

^±(^) = y"^ — ^S(uj-e- AsTnujp) (163) 

■^-^ nl 

n=0 

where -|- refers to absorption spectrum and — to emission spectrum, a = ^ 9q/^p ^^'^ ^^ = ^^p 
is the shift in core energy due to the interaction with the plasmon field. It is assumed that the 
plasmon excitations have no dispersion although this assumption is not necessary. The spectra 
consist therefore of the main quasiparticle peak at lu — e + As and a series of plasmon excitations 
at multiples of the plasmon energy below the quasiparticle peak which is in accordance with 
experiment. This is in contrast to the GH^ spectra which only has one plasmon excitation located 
at too high energy, typically 1.5 tOp below the quasiparticle peak (Hedin, Lundqvist, and Lundqvist 
1970). More recently the cumulant expansion method was applied to a model Hamiltonian with 
electron-boson interaction and the cumulant was calculated to higher order (Gunnarsson, Meden, 
and Schonhammer 1994). 

In the cumulant expansion approach, the Green function for the hole (t < 0) is written as (Lan- 
greth 1970, Bergersen, Kus, and Blomberg 1973, Hedin 1980) 

G (fc, t < 0) = i9 i-t) {N\cl{0)ck{t)\N) 

= ^6l(-t)e-*^''•*+^'('=^■*) (164) 

and the hole spectral function is 

A{k,Lu < fi) = -ImG (k, uj < ^) 

= - J ^dte^-\N\cl{0)ck{t)\N) 

= -Imt f die*'^*e-^^^*+<^'('='*) (165) 

T^ J-oo 

where k denotes all possible quantum labels and C^ {k, t) is defined to be the cumulant. Expanding 
the exponential in powers of the cumulant we get 



G{k,t)^Go{k,t) 



l + G^fc,t) + i[G^fc,f)]2 + 



(166) 



where Gq {k, t) — i exp {—iskt) . In terms of the self-energy, the Green function for the hole can be 
expanded as 

64 



G = Go + GoSGo + GoSGoSGo + . . . (167) 

To lowest order in the screened interaction W, the cumulant is obtained by equating 

GoG'' = GoSGo (168) 

where S — T,cw = iGqW. If Go corresponds to, e.g. Glda, then E = Y^gw — V'™. The first-order 
cumulant is therefore (Hedin 1980, Almbladh and Hedin 1983, Aryasetiawan, Hedin, and Karlsson 
1996) 



/OO /'OO 

dt' / dre'^'^J^ik,' 



(169) 



The cumulant may be conveniently divided into a quasiparticle part and a satellite part: G'* 
Cqp + Cs where 



with 



C^p{k, t) = {iak + 7fe) + {-iAskt + f]k)t 

/M „i(efc-w-i(5)t 
dc^- -2r(fc,^) 
-OO [Ek — U! — 10) 

ai](fc,tj) 



Jafc + Ik = 



duj 



(170) 
(171) 

(172) 



Asfc = P 



-p / ; \ 

dw ^ — , ryfc = 7rr(fc,efe) = |ImI](k,ek)| 



(173) 



T^fc is the inverse life-time of the quasiparticle and r(fc, oj) is the spectral function of the self-energy 
which is proportional to ImS(k, a;). A similar derivation can be carried out for the particle Green 
function 



The result is 

CP {k,t) = ~i/\ekt - f^kt + 



dJ:{k,Lu 



duj 



duj -T (fc, Lo) 



(efc - cj + id) 



It is physically appealing to extract the quasiparticle part from the Green function: 
G'^^P (k, t) = iB i-t) e^"^+7.e(-*^'=+'"=)* , Ek = ek + Ask 



(174) 



(175) 



(176) 



The spectral function for this quasiparticle can be calculated analytically (Almbladh and Hedin 
1983): 



Aqp {k,Lu < ^) 



e '"' r]k cos ttk ~ {lo - Ek) sin ak 
{iv-Ekf + vl 



(177) 



From equation ( |l65|) we have {N\cl{0)ck{t)\N) = Q- iekt+c"-{k,t) j-q^. t < q. By analytically contin- 
uing to f > the spectral function in equation (|165|) can be rewritten as 



(178) 



1 r°° 
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where for i > we have C'^{k, t) = C^*{k, —t). The total spectra can be written as a sum of Aqp 
and a convolution between the quasiparticle and the satellite part: 



A{k,Lu) = AQp{k,Lu) + 



2tt 



^^^i<^t^(-iEk+Vk+C'\kfi))t 



,CS(fc,t) _ I 



= Aqp (fc, uj) + Aqp {k, u) * As (fc, uj) 



(179) 



where 



As{k,u;) =^ /'dte'-*{e^S(M)_iJ 



(180) 



The second term Aqp * As is responsible for the satellite structure. The Fourier transform of Cg 
can be done analytically (Aryasetiawan, Hedin, and Karlsson 1996) 



C^(fc,tj<0) = 



r{k,ek+u;)-r{k,Sk)-u;r{k,Sk) 



(181) 



As follows from equations (177) and ( |l79| ), the quasiparticle energy in the cumulant expansion 
is essentially determined by Ek, which is the quasiparticle energy in the GWA. 



Cjtgh- = 



-» '-—»' 1 f 1— -» 1 ►- 



VJcE/A; 



-»• — i— -» — I f — I— -» — I ►- 



, r^ . 



FIG. 20. Diagrammatic expansion for the Green function 
to second order in the GWA and the cumulant expansion 
respectively. The solid lines represent non-interacting Green 
functions Go, and the wiggly lines represent the screened 
interaction W . 
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FIG. 21. (a) The experimental spectral function of Na (dots). The solid line is a synthetic spectrum 
obtained by convoluting the density of states from a band structure calculation and the experimental 
core level spectrum. BG is the estimated background contribution. From Steiner, Hochst, and Hiifner 
(1979). (b) The calculated total spectral function of Na for the occupied states. The solid and dashed 
lines correspond to the cumulant expansion and the GWA, respectively. After Aryasetiawan, Hedin, and 
Karlsson (1996). 
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By comparing the diagrammatic expansions in the GWA and the cumulant expansion we can 



get some idea about the vertex corrections. In figure 20 the Green function diagrams are shown 
to second order in the screened interaction, which should be sufficient for our purpose. The 
cumulant expansion diagrams are obtained by considering the three possible time-orderings of the 



integration time variables t' in C^ {k,t) with C^ {k,t) given by equation (169). The cumulant 
expansion contains second-order diagrams which are not included in the GWA. It is these additional 
diagrams that give rise to the second plasmon satellite and they are quite distinct from the second- 
order diagram common to both approximations. The interpretation of the latter diagram is that a 
hole emits a plasmon which is reabsorbed at a later time and the hole returns to its original state 
before plasmon emission. This process is repeated once at a later time. Thus there is only one 
plasmon coupled to the hole at one time. In contrast, the other two diagrams, not contained in 
the GWA, describe an additional plasmon emission before the first one is reabsorbed, giving two 
plasmons coupled to the hole simultaneously. Similar consideration can be extended to the higher- 
order diagrams. If self-consistency is taken into account then the second second-order diagram is 
also included in the GWA. 

The cumulant expansion contains only boson-type diagrams describing emission and reabsorb- 
tion of plasmons but it does not contain diagrams corresponding to interaction between a hole and 
particle-hole pairs. For this reason, the cumulant expansion primarily corrects the satellite descrip- 
tion whereas the quasiparticle energies are to a large extent determined by the GWA as mentioned 
before. Interaction between hole-hole and particle-hole is described by the ladder diagrams which 
in a Hubbard model study were found to improve the low-energy satellite (Verdozzi, Godby, and 
Holloway 1995). 

The cumulant expansion was applied recently to calculate the photoemission spectra in Na and 
Al (Aryasetiawan, Hedin, and Karlsson 1996). The experimental spectra consist of a quasiparticle 
peak with a set of plasmon satellites separated from the quasiparticle by multiples of the plasmon 
energy. The spectra in the GWA shows only one plasmon satellite located at a too high energy, 
approximately 1.5 LUp below the quasiparticle which is similar to the core electron case. The 
cumulant expansion method remedies this problem and yields spectra in good agreement with 
experiment regarding the position of the satellites. The relative intensities of the satellites with 
respect to that of the quasiparticle are still in discrepancy. This is likely due to extrinsic effects 
corresponding to the interaction of the photoemitted electron with the bulk and the surface on 
its way out of the solid resulting in energy loss. These are not taken into account in the sudden 
approximat ion . 

When applied to valence electrons with band dispersions the cumulant expansion does not yield 
the exact result anymore as in the core electron case. Surprisingly, the numerical results show 
that the cumulant expansion works well even in Al with a band width of ~11 eV. Considering its 
simplicity, it is a promising approach for describing plasmon satellites. 



VIII. SUMMARY AND CONCLUSIONS 

The GWA has been applied by now to a large number of systems, ranging from atoms, simple 
metals, semiconductors, transition metals, clusters, and surfaces and interfaces. In practically all 
of these systems, the GWA improves the quasiparticle energies relative to the LDA eigenvalues. 
The reason for the success of the GWA may be understood qualitatively by the fact that it is 
correct in some limiting cases as described in the introduction. The GWA includes an important 
physical ingredient in extended systems, namely screening or polarization of the medium which 
is absent in the HFA. It is well known that the neglect of screening leads to unphysical results 
in metals such as zero density of states at the Fermi level and in semiconductors and insulators 
to too large band gaps. Even in atoms, the inclusion of polarization effects lead to a significant 
improvement on the HF eigenvalues. Since screening is a common feature in all electronic systems, 
it is perhaps not surprising that the GWA works in a wide variety of materials. It is also known 
empirically that second order perturbation theory often takes into account most of the physical 
effects, in particular the shifting of quasiparticle energies. 
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Despite of the success of the GWA, it has naturally some shortcomings. One of these is related 
to satellite structure in the photoemission spectra. Since the GWA describes a coupling of the 
electrons to one plasmon excitation, represented by W, multiple plasmon excitations observed in the 
alkalis are clearly beyond scope of the GWA. This problem is remedied by the cumulant expansion 
theory described in section VII-C. Apart from plasmon-related satellites, there are also satellite 
structure originating from short-range interactions. This type of satellites appears in strongly 
correlated materials containing 3d or 4f orbitals. The GWA is based on the RPA screening which 
takes into account the dominant part of long-range (small momentum) screening. Short-range 
(large momentum) intrasite interactions of multiple holes usually present in strongly correlated 
systems are therefore not well described by the GWA. Here a T-matrix type approach may be 
appropriate as has been shown by model calculations. Apart from problems with satellite, there 
are also in some cases discrepancies in the quasiparticle energies. For example, the band width in 
Na is « 10 % too large within the GWA compared with experiment but these discrepancies are 
relatively small. 

Another shortcoming is the absence of spin-dependence in the screened interaction W since 
the screening is purely Coulombic. The spin dependence enters only through the Green function. 
One would therefore expect some problems when applying the GWA to magnetic systems where 
spin-spin correlations are important as indicated by GW calculations on transition metal atoms. 
This area of research has not been explored extensively and it would require inclusion of vertex 
corrections to take into account spin-spin correlations. A T-matrix approach may be a first step 
in this direction. 

Most GW calculations performed so far are not self-consistent, i.e. the Green function used to 
calculate the self-energy is not equal to the Green function obtained from the Dyson equation with 
the very same self-energy. Only very recently such self-consistent calculations were performed for 
the electron gas. The results turn out to be worst than the straight GW calculations (one iteration) 
and clearly show a cancellation between the effects of self-consistency and vertex corrections. In 
some way it is a blessing since fully self-consistent calculations are numerically difficult to perform 
for real systems. One interesting aspect of the self-consistent calculations is that the total energies 
are in very good agreement with the QMC results. A recent finding shows that when the total 
energies are calculated within the Luttinger-Ward functional and its extension (Almbladh, von 
Barth, and van Leeuwen 1997) using approximate G and W the results are almost equally good, 
which circumvents the need for self-consistent calculations (Hindgren 1997). This could be due to 
the variational property of the functional and the fact that the self-consistent GWA is conserving 
in the sense of Baym and Kadanoff. Calculating total energies within the GWA could be an 
alternative to the QMC method, particularly for extended systems. 

Due to computational difficulties, applications of the GWA to large and complex systems are 
still not feasible. Simplified GW schemes, which are computationally efficient and yet maintain the 
accuracy of the full calculations, are therefore very desirable. Many schemes have been proposed 
but most of them are designed for semiconductors. While they give reliable band gaps, details 
of the band structure are not fully accounted for. Reliable schemes must probably take into 
account non-locality as well as energy dependence of the self-energy. With efficient schemes, many 
interesting problems can then be tackled. These include chemisorption at surfaces, 3d impurities in 
semiconductors, interfaces, band off-sets in heterojunctions, and exotic materials such as fuUerenes 
in their diverse forms. 
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